% code to simulate diffusion process with   Dirichlet BC's
%via Method of lines
%
function Diff_sim
 global alpha sc scu scl bc

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
 alpha = D/dx^2;  %The time constant
 
 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
 V0 = exp(-(X-dx*N/2).^2./(2*dx));  % initial data is a gaussian
 
tstep = 0.01; 
t_end = 0.5; %total time to run simulation
 
%specify the output points
tspan = [0:tstep:t_end];

 
% integrate the system of ode's:
% pick an ode integrator from the various Matlab options
[T,S] = ode45(@deRHS,tspan, V0, odeset('maxstep',1));  

% plot the solution at each time step
figure(1)
for j = 1:length(T)
    plot([0;X;L],[u0,S(j,:),uL],'*')
    ylabel('u(X,t)','fontsize',20)
    xlabel('X','fontsize',20)
    axis([0 L 0 1])
    pause(.1)
end
figure(2)
plot(T,S(:,N/2),'linewidth',2)
xlabel('t','fontsize',20)
ylabel('u(N/2,t), u(N/2+1,t) ','fontsize',20)

 
%the right hand side for ode simulation:
function s_prime=deRHS(t,u)
global alpha sc scu scl bc

 %the spatial difference operator:
 
Fu = alpha*(-sc.*u+scl.*[0;u(1:end-1)]+scu.*[u(2:end);0] +bc) ;
s_prime = Fu;

