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

 L = 100; %length of the domain
D = 1; % the diffusion coefficient
 
 N = 80; %number of interior grid points  (= number of unknowns u_j)
 dx = L/(N); % grid size 
  
 sc =  2*ones(N+1,1) + [0;zeros(N-1,1);0]; % 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
  
 alpha = D/dx^2;  %The time constant
 h=dx;
  
 X = [0:N]'*dx;  %  grid points
 V0 = exp(-(X-dx*N/2).^2./(2*dx));  % initial data is a gaussian
 V0=exp(-0.05*X.^2);
tstep = 0.5; 
t_end = 25.; %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] = ode23(@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
 
% estimate the speed
contour(T,X,S',[0.5 0.5],'linewidth',2)  % make a contour plot of the V = 0.5 level surface
 ylabel('x', 'fontsize',20)
  xlabel('t', 'fontsize',20)
  
  % now a little post processing of S to find the speed of travelling
  % waves
  M = length(T);
   xf = zeros(M,1);
  for j = 1:M
      if (max(S(j,:))>=0.5)
      xj = max(find(S(j,:)>=0.5))
      % estimate the x value where S = 0.5
      if(xj<N)
   
      Xf(j) = X(xj) +(0.5-S(j,xj))*h/(S(j,xj+1 )-S(j,xj ));  %This estimates where V=0.5
      end
      end
      
  end
  
  jmin = min(find(Xf>0));
  tmin= T(jmin);
  xmin = Xf(jmin);
  
  jmx = max(find(Xf<L));
  tmx= T(jmx);
  xmx = Xf(jmx);
  
  rngcut = 15;
  rnge = [jmin+rngcut:jmx-rngcut];  % shorten the range a bit to get a better straight line
  
  p = polyfit(T(rnge),Xf(rnge)',1)
  p(1) % is an estimate of the speed
  figure(4)
  plot(T(rnge) ,Xf(rnge) ,T(rnge),p(1)*T(rnge)+p(2),'--','linewidth',2)
  
  
 
%the right hand side for ode simulation:
function s_prime=deRHS(t,u)
global alpha sc scu scl   

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

