 % This file solves the diffusion equation with Dirichlet BC's
 %using the backward Euler method
clear
 
L = 1; %length of the domain
D = 1; % the 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

 
 sc =  2*ones(N,1) ; % diagonal of difference operator
 scu = [ones(N-1,1)]; %upper diagonal for difference operator 
 scl = [ones(N-1,1)]; %lower diagonal for difference operator
  
M = 20;  % number of time steps
%pick a time step:
dt = dx/2; % time step  %this is much bigger than for forward euler!
scal = D*dt/dx^2;

 
 V = exp(-(X-dx*N/2).^2./(2*dx));  % initial data is gaussian
 
kstep = 1;
t=0;

% uses backward Euler to solve the diffusion equation

% set up the implicit matrix
 
Am=diag(1+scal*sc)-scal*diag(scu,1)-scal*diag(scl,-1);

% The algorithm
 
for n=2:M
for k = 1:kstep
 
V = Am\(V+scal*bc);
 
t=t + dt;
end
figure(1)
  plot([0;X;L],[u0;V;uL],'linewidth',2 )
  axis([0 L  0 1])
 xlabel('x','fontsize',20)
 ylabel('u','fontsize',20)
   pause(0.05)

end
  
