%code1 for lecture 2; need u_0(x) and f(x, t) dt=input(' dt= '); T=input( 'T= '); alpha=input( 'alpha= '); for N=2:12 d1=DM1(N); d2=d1*d1; A=d2(2:N,2:N); %initial condition j=[1:1:N-1]; x=[cos(j*pi/N)]'; u=sin(pi*x); time=0.0; while time <= T G=A*u; u=u+alpha*dt*G; G=(-1+2*alpha-2*alpha^2)*G+A*u; u=u+dt*G/(2*alpha); time=time+dt; end %exact solution exact=exp(-pi^2*time)*sin(pi*x); error=abs(exact-u); xx(N-1)=N; err(N-1)=max(error); ss=max(error) end fprintf(1, '%9.0f %10.2e \n', [xx; err]);