% code to simulate diffusion equation via Method of lines
%with  Neumann or Robin BC's
function Diff_sim
 global alpha sc scu scl bc

 L = 10; %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
 

 alpha = D/dx^2;  %The time constant
 
  
 X = [0:N]'*dx  %  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(X, S(j,:) ,'*')
    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;

