%An example of a matrix D for which is is hard to compute the inverse but easy to solve equations D x = b

clear;
% define matrix

% size of the matrix
n=40;

% off-diagonal part
B=zeros(n,n);
for i=1:n-1
  B(i,i+1)=1;
end

% diagonal part
lam=-0.1;
D=B-lam*eye(n);


% define first and last unit vector
en=[zeros(1,n-1) 1]';
e0=[1 zeros(1,n-1)]';

% compute inverse and measure error
Dinv=inv(D);
errorinverse=max(max(Dinv*D-eye(n)))

% solve equation D u = D(:n) 
% which has the (obvious) solution en
% using

% ...the previously computed inverse
errorinvnthcolumn=max(Dinv*D(:,n)-en)

% ...gauss elimination
errorlrnthcolumn=max(D\D(:,n)-en)


% ...an LU decomposition
[L,U]=lu(D);
errorLUinv=max((U\(L\D(:,n)))-en)
% note that we do not compute inv(U)!
