% code to simulate bistable eqn via Method of lines

%  
function Diff_sim
 global alpha sc dg
 
 set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 

cmykflg=0;
%if cmykflg=1, this code will make cmyk figure files
%cmykflg = 1;

% first simulation - propagation failure
 
 dglist=[0.018,0.02];
 tlist=[300,800];
 
 for jd = 1:length(dglist)
 dg = dglist(jd)
 t_end=tlist(jd);
 
alpha = 0.25 ; % the threshold
 N = 14; %number of cells    
 sc = [1;2*ones(N-2,1);1];
 u0 = zeros(N,1);
 u0(1:5) = 0.5;  %the first five cells are above threshold
 
tstep = 0.5;
 
%specify the output points
tspan = [0:tstep:t_end];

[T,S] = ode23(@deRHS,tspan, u0, odeset('maxstep',1));  

figure(jd)
for j = 1:N-1
    plot(T,S(:,j),'linewidth',2)
     
   hold on
end
 hold off
 
xlabel('t','fontsize',20)
ylabel('u_n(t)','fontsize',20)
axis([0 t_end 0 1])
switch jd
    
    case 1
        
        
text(25,0.13,'n=0','fontsize',20)
text(25,0.94,'n>0','fontsize',20)
text(25,0.04,'n<0','fontsize',20)
  if(cmykflg==1)
  print('../../figs_c/chapt_8/discrete_fail','-deps','-cmyk')
   end 
 figure(length(dglist)+1)
 shift = 6;
 plot([1-shift:N-shift],S(end,:),'*','linewidth',2)
 hold on
 plot([1-shift:N-shift],S(end,:),'--')
 hold off
 axis([1-shift N-shift 0 1])
xlabel('n','fontsize',20)
ylabel('u_n(\infty)','fontsize',20)
 if(cmykflg==1)
  print('../../figs_c/chapt_8/discrete_inf_prf','-deps','-cmyk')
   end 
    case 2
  if(cmykflg==1)
  print('../../figs_c/chapt_8/discrete_prop','-deps','-cmyk')
   end       
end
% 
% % Now find the speed:
[T,n] = meshgrid(T,1:N);
figure(4)
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(5)
plot(spx,spy,spx,spest,'--')
 xlabel('\tau', 'fontsize',20)
  ylabel('n', 'fontsize',20)
  
 end
 
%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;

