  clear all; clc
  D=1; mu=1; A=1; B=1; 
  nx=21; xl=1; dx=0.05; x=[0:0.05:xl]; nx2=11;
  nt=50; nout=10; h=0.001; t=0;
  for i=1:nx; w(i)=A*exp(1/D*(mu^2*t-mu*x(i)))+B; end
  for i1=1:nout
    for i2=1:nt
      for i=1:nx
        if(i==1)w(1)=A*exp(1/D*(mu^2*t-mu*x(1)))+B; wt(1)=0;   
        elseif(i==nx)wf=w(nx-1)+(2*dx)*A*(-mu/D)*exp(1/D*(mu^2*t-mu*x(nx))); 
                     wt(nx)=D*(wf-2*w(nx)+w(nx-1))/dx^2; 
        else wt(i)=D*(w(i+1)-2*w(i)+w(i-1))/dx^2; end 
      end 
      w=w+wt*h; t=t+h;
    end  
  fprintf('t = %7.4f   w(x=xl/2,t) = %7.4f   wa(x=xl/2,t) = %7.4f\n',...
          t,w(nx2),A*exp(1/D*(mu^2*t-mu*x(nx2)))+B);
  end
      
