 % This file solves the diffusion equation using Forward Euler algorithm
clear
 
L = 1; %length of the domain
D = 1; % the diffusion coefficient
del0 =  0; % porosity or reaction rate at x=0 (for Robin condition)
delL = 0; % porosity or reaction rate at x=L (for Robin condition)
% Robin BC at x=0:  -Du_x + del0 u0 = U0 at x=0
%Robin BC at x=L: Du_x +delL u_N = UL  
U0 =  0;
UL = 0;

 N = 20; %number of interior grid points  (= number of unknowns u_j)
 dx = L/(N ); % grid size 
 

 rsc0 = 2*del0*dx/D; %the robin condition diagonal factor 
 rscL = 2*delL*dx/D;  %the robin condition diagonal factor
 sc =  2*ones(N+1,1) + [rsc0;zeros(N-1,1);rscL]; % diagonal of difference operator
 scu = [2;ones(N-1,1);0]; %upper diagonal for difference operator 
 scl = [0;ones(N-1,1);2]; %lower diagonal for difference operator
 bc = 2*dx/D*[U0;zeros(N-1,1);UL]; % for inhomogeneous BC's
 
 
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;

X = dx*(0:N)'; % the N+1 grid points
 % specify initial data
V = exp(-(X-dx*N/2).^2./(2*dx));  %initial data is gaussian
figure(1) % plot initial data
  plot(X,V )
  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
 % here is the discretization of the laplacian:
V_t =  scal * (scu.*[V(2:N+1 );0] - sc.*V + scl.*[0;V(1:N )]+bc); 
 
% update V
V = V+V_t;
 
t=t + dt;
end
figure(1)
  plot(X,V,'linewidth',2 )
  axis([0 L  0 1])

   pause(0.05)

end
  
