 % This file solves the diffusion equation with Dirichlet BC's 
 %using Forward Euler algorithm
clear
 
L=1; %length of domain
D = 1; % diffusion coefficient
u0 = 0; % boundary value at x=0
uL = 0 ; % boundary value at x=L
 N = 20; %number of interior grid points  (= number of unknowns u_j)
 dx = L/(N+1); % grid size
 
 bc = [u0;zeros(N-2,1);uL]; % boundary contribution to difference operator
 sc =  2*ones(N,1) ; % diagonal of difference operator
 scu = [ones(N-1,1);0]; %upper diagonal for difference operator 
 scl = [0;ones(N-1,1)]; %lower diagonal for difference operator
 
 X = [1:N]'*dx;  % interior grid points

 
 
M = 25;  % number of time steps
%pick a time step:
dt =  dx^2/(4*D); % time step  % must have dt<h^2/(2D)
 
scal = D*dt/dx^2;

 % specify initial data
V = exp(-(X-dx*N/2).^2./(2*dx));  %initial data is gaussian
figure(1) % plot initial data
  plot([0;X;L],[u0;V;uL],'linewidth',2 )
  axis([0 L  0 1])
 
kstep = 10;
t=0;

% uses Forward Euler to solve the diffusion equation

% The algorithm

for n=2:M
for k = 1:kstep
 
V_t =  scal * (scu.*[V(2:N );0] - sc.*V + scl.*[0;V(1:N-1 )]+bc)  ;% with Neumann BC's
 
% update V
V = V+V_t;
 
t=t + dt;
end
figure(1)
  plot([0;X;L],[u0;V;uL],'linewidth',2 )
  axis([0 L  0 1])

   pause(0.05)

end
  
