  clear all; clc
  c=1; 
  nx=201; xl=50; dx=0.5; x=[-xl:dx:xl]; 
  nt=4000; nout=3; h=0.0025; t=0;
  for i=1:nx; w(i)=exp(-0.05*x(i)^2); wt(i)=0;...
              wplot(1,i)=w(i); waplot(1,i)=exp(-0.05*x(i)^2); end 
  for i1=1:nout
    for i2=1:nt
      for i=1:nx
        if(i==1)     w( 1)=0; wt( 1)=0; wtt( 1)=0;   
        elseif(i==nx)w(nx)=0; wt(nx)=0; wtt(nx)=0; 
        else wtt(i)=c^2*(w(i+1)-2*w(i)+w(i-1))/dx^2; end 
      end 
      w=w+wt*h; wt=wt+wtt*h; t=t+h;
    end 
  for i=1:nx 
     wplot(i1+1,i)=w(i);  
    waplot(i1+1,i)=0.5*(exp(-0.05*(x(i)-t)^2)+exp(-0.05*(x(i)+t)^2));
    fprintf('\n %6.2f %6.2f,%8.4f %8.4f',t, x(i),wplot(i1+1,i),waplot(i1+1,i));
  end    
  end
  plot(x,wplot,'-',x,waplot,'o'); xlabel('x'); ylabel('w(x,t)'); 
  title('w(x,t), t=0,10,20,30, solid - num, o - anal')
      
