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(