%Explicit numerical method for u_t = k u_xx

clear;  %  gets rid of old data

L = 50;  % domain size
eps = 2.0;   % "diffusion" constant

N = 200; % number of spatial mesh points
dx = L/(N+1);

dt = 0.4*dx^2/eps^2;  %  timestep is limited by  1/2 (dx/eps)^2; Why?

%setup grid
x = [dx: dx: L-dx]';   % x is the column vector x_1, ..., x_N-1
                       % the ' denotes a matrix transpose

% setup initial data:  u is a N-1 size column vector which
% will always contain the solution at the present time step

left = 2;    %left and right boundary conditions
right = 0;

%u =  (rand(N,1)).*(x > 10).*(x < 15);       % random pulse initial data
u = 1 - tanh( (x - 10)/eps);  % this initial data looks like a smooth step between 2 and 0, centered at x=10

% now step forward
time=0;
while (time < 100),
   % these commands will plot the solution and pause
   % briefly before proceeding to the next time-step.
   % to make things go slower, increase the pause time
   plot(x,u);
   axis([0 50 -1 3]);
   text(40,2.5,['time=' num2str(time)]);
   pause(.05);

   % do points next to boundary separately
   unext(1) =  u(1) + dt*( eps^2/dx^2*(-2*u(1) + u(2) + left )  - u(1)*(u(2)-left)/(2*dx)  );
   unext(N) = u(N) + dt*( eps^2/dx^2*(-2*u(N) + u(N-1) + right)   - u(N)*(right - u(N-1))/(2*dx) );
   %then do rest (note the "vectorization")
   unext(2:N-1) = u(2:N-1) + dt*(eps^2*(u(3:N)+u(1:N-2)-2*u(2:N-1))/dx^2 +... 
                    - u(2:N-1).*(u(3:N)-u(1:N-2))/(2*dx) );
   
   u = unext;          % now u is equal to the updated solution
   time = time + dt;
end;


