clear all;

M=[16 32 64 128 256 512];
e=zeros(1,length(M));
a=1;
for k=1:length(M)
    N=M(k);
    h=1/N;
    x=[h:h:1-h];
    A =  diag((-2/h^2-a)*ones(N-1,1))+diag( ones(N-2,1)/h^2,1) ...
        + diag(ones(N-2,1)/h^2,-1); 
    [B, D, C]=crout(A);
    b=sin(pi*x);
    %solve LY=b, where Y = U*(uh)
    Y=zeros(N-1,1);
    Y(1) = b(1);
    for i = 2:N-1
        Y(i) = b(i)-B(i)*Y(i-1) ;  % B(i) is L(i,i-1), i = 2, ..., N-1 
    end
    %solve U*(uh) = Y
    uh=zeros(N-1,1);
    uh(N-1) = Y(N-1)/D(N-1);       % D(i) is U(i,i), i = 1, ..., N-1
    for ii = N-2:-1:1
        uh(ii) = ( Y(ii)-C(ii)*uh(ii+1) )/D(ii); % C(i) is U(i,i+1)
    end
    ue=-sin(pi*x)/(pi^2+a);
    e(k)=max(abs(ue-uh'));
end
figure(1)
axis normal,
loglog(1./M,e,'x')
title('HW10, problem 4(c)')
xlabel(' log(h)')
ylabel('log( max|uh-ue| )')