%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]);
