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*N^2-a)*ones(N-1,1))+diag((N^2)*ones(N-2,1),1)+diag((N^2)*ones(N-2,1),-1); [L ,U]=al67(A); b=sin(pi*x); %solve LY=b, where Y = U*(uh) Y=zeros(N-1,1); for i = 1:N-1 tmp = 0; for j =1:i-1 tmp = tmp + L(i, j)*Y(j); end Y(i) = (b(i)-tmp)/L(i, i); end %solve U*(uh) = Y uh=zeros(N-1,1); for i = 1:N-1 tmp = 0; for j =1:i-1 tmp = tmp + U(N-i, N-j)*uh(N-j); end uh(N-i) = (Y(N-i)-tmp)/U(N-i, N-i); 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| )')