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); [L ,U]=al67(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)-L(i,i-1)*Y(i-1) )/L(i, i); end %solve U*(uh) = Y uh=zeros(N-1,1); uh(N-1) = Y(N-1)/U(N-1,N-1); for ii = N-2:-1:1 uh(ii) = ( Y(ii)-U(ii,ii+1)*uh(ii+1) )/U(ii, ii); 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| )')