N=20; A = diag(5.1*ones(N,1))+diag(2*ones(N-1,1),1) + diag(2*ones(N,1),-1); [B, D, C]=crout(A); bb=ones(20,1); uh=zeros(N,1); z(N) = bb(N)/D(N); for ii = N-1:-1:1 z(ii) = ( bb(ii)-C(ii)*z(ii+1) )/D(ii); % C(i) is U(i,i+1) end max(abs(