%Explicit numerical method for u_t = k u_xx

clear;  %  gets rid of old data

L = 10;  % domain size
k = 3;   % diffusion constant

dt = .001; % time step size
M = 300; % total number of time steps

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

s = dt*k/dx^2;  %derived parameter that appears in the numerical method

%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

u = x.*(10-x);  % initial data u = x*(10-x). 


% now step forward
time=0;
for j = 1:M,
   % 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 10 -10 25]);
   text(0.1,20,['time=' num2str(time)]);
   pause(.05);

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