h=0.1; k=0.0005; gamma=k/h^2; m=1/h; A=zeros(m-1,m-1); A(1,1)=1-2*gamma; x=linspace(h,1-h,m-1)'; for i=2:m-1, A(i,i)=1-2*gamma; A(i-1,i)=gamma; A(i,i-1)=gamma; end w=sin(pi*linspace(h,1-h,m-1))'; for i=1:(0.5/k), w=A*w; true=exp(-pi^2*i*k)*sin(pi*x); plot(x,true,x,w); pause(.005); end