function A = Aoneseps(n)
% A = Aoneseps(n);
%
% Calculates a matrix A of size nxn for which Gaussian
% elimination is unstable; and, results with different algorithms are
% displayed.

disp(' ')
disp(' ')
disp('starting new random example')
disp(['first generate special A and random x, and b=Ax, for n= ',...           
                  num2str(n)]);
A=-tril(ones(n));
A(1:n+1:end)=ones(n,1);
A(:,end)=ones(n,1);
xhat=randn(n,1);
b=A*xhat;

disp(' ')
disp('Solve unperturbed problem using x = A\b:');
pause(.1)

tic, x = A\b; toc
res = norm(A*x-b);  disp(['residual  norm(A*x-b) ',num2str(res)]);
err = norm(x-xhat); disp(['error  norm(x-xhat) ',num2str(err)]);



E=randn(n);
Aeps=A+   100*eps*  (  E/ (max(max(abs(E))))  );
disp('Now *****add small random perturbation***** to A --> Aeps')
disp('Solve by xeps = Aeps\b:');
pause(.1)

tic, xeps = Aeps\b; toc
res = norm(Aeps*xeps-b);  disp(['residual norm(Aeps*xeps-b) ',num2str(res)]);
res = norm(A*xeps-b);  disp(['residual  norm(A*xeps-b) ',num2str(res)]);
err = norm(xeps-xhat); disp(['error  norm(xeps-xhat) ',num2str(err)]);
