% code to simulate bistable eqn via Method of lines

%  
% code to find speed for discretely coupled cellsfunction Diff_sim
function discrete_speed
global alpha sc dg
 
 set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 
 % set up the parameters
alpha = 0.25 ; % the threshold for the bistable cubic function

N = 100; %number of cells  ; this affects the accuracy of the calculation  
 sc = [1;2*ones(N-2,1);1];  % this is the spatial discretization operator for MOL
 
 %dg is the coupling coefficient assumed to be D/dx^2;
 %Take D=1;
 % for dg <0.02, propagation fails for this value of alpha
 Nj = 50;  % number of points
 dx = (7.07-0.5)/Nj;
 for jg = 1:50
     DX = 0.5+ (jg-1)*dx;
     dg = 1/DX^2;
 u0 = zeros(N,1);
 u0(1:5) = 0.9;  %the first five cells are above threshold
 
tstep = 0.5;
t_end =500;

%specify the output points
tspan = [0:tstep:t_end];

[T,S] = ode23(@deRHS,tspan, u0, odeset('maxstep',1));  
 
figure(1)
for j = 1:N-1
    plot(T,S(:,j),'linewidth',2)
   hold on
end
 
xlabel('t','fontsize',20)
ylabel('u_n(t)','fontsize',20)

hold off
% 
% % Now find the speed:
figure(2)
[T,n] = meshgrid(T,1:N);
sfc=contour(T,n,S',[0.5 0.5],'linewidth',2);
sfcl=length(sfc(1,:));
% to get a good estimate of the speed, this curve may need to be truncated
% some
spx = sfc(1, 2:end);
spy = sfc(2, 2:end);
 
  % now use sfc to estimate the speed
  % this must be done interactively
  
  p=polyfit(spx,spy,1) %p(1) is the estimated speed
  spest = p(2)+p(1)*spx;
  
  figure(3)
plot(spx,spy,spx,spest,'--')
hold on

 xlabel('\tau', 'fontsize',20)
  ylabel('n', 'fontsize',20)
  dj(jg) = dg;
  spd(jg) =  p(1);
 end
 hold off
 
 c0=sqrt(2)*(0.5-alpha);
 figure(14)
 plot(dj, spd,dj, sqrt(dj)*c0,'--')
 xlabel('d','fontsize',20)
 ylabel('Speed = number/time','fontsize',20)
 
 figure(15)
 % set d = D/dx^2, with D=1;
 %
 cc = sqrt(2)*(0.5-alpha);  % continuous speed
 
 dx = sqrt(1./dj);
 plot(dx,dx.*spd,dx, cc*ones(1,length(dx)),'--','linewidth',2)
 xlabel('\Delta x','fontsize',20)
 ylabel('Speed = Distance/time','fontsize',20)
 
%the right hand side for ode simulation:
function s_prime=deRHS(t,s)
global alpha sc dg
 
      Fu =  dg*(-sc.*s+[0;s(1:end-1)]+[s(2:end);0]) +...
          s.*(s-alpha).*(1-s);;
s_prime = Fu;

