% code to simulate diffusion process via Method of lines with no flux BC's
function Diff_sim
 global alpha sc
 
 set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);

alpha = 0.1;  %The time constant
  
 N = 52; %number of cells (needs to be an even number)    
 sc = [1;2*ones(N-2,1);1];
 u0 = zeros(N,1); %initial data
 u0(N/2) = 1; 
 
tstep = 2; 
t_end = 100; %total time to run simulation
 
%specify the output points
tspan = [0:tstep:t_end];

 
% integrate the system of ode's:
[T,S] = ode23s(@deRHS,tspan, u0, odeset('maxstep',1));  

% plot the solution at each time step
figure(1)
for j = 1:length(T)
    plot(S(j,:),'*')
    ylabel('u(n,t)','fontsize',20)
    xlabel('n','fontsize',20)
    axis([0 N 0 1])
    pause(.1)
end
figure(2)
plot(T,S(:,N/2),T,S(:,N/2+1),'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,s)
global alpha sc
 
Fu = alpha*(-sc.*s+[0;s(1:end-1)]+[s(2:end);0]);  % this is the discretized diffusion operator
s_prime = Fu;

