% code to simulate Gray-Scott diffusion -reaction process via Method of lines
function Diff_sim
 global N sc   rho kv L du dv
 
 set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 

 

 %first make a plot of parameter region:
 rho0 = [0:0.001:0.25];
k0 = sqrt(rho0)/2-rho0;

  % set the parameters
rho = 1/16;
 kv = 0.02;
 
 figure(10) % not used in the text, but is used for exercise 12.2
 
 plot(rho0,k0,rho,kv,'*','linewidth',2)
 xlabel('\rho','fontsize',20)
 ylabel('k','fontsize',20)
 text(0.05,0.03,'bistable','fontsize',20)
 text(0.15,0.055,'monostable','fontsize',20)
  
%Make a phase plane portrait for the ode dynamics
u = [0.001:0.001:1];
% null clines
v1 = sqrt(rho*(1-u)./u);
v2 = (rho+kv)./u;

% integrate the ode
tstep = 1;
t_end = 50;
 
%specify the output points
tspan = [0:tstep:t_end];

s0 = [1;0.095];  %initial data for the ode
[T,S] = ode23(@deRHS,tspan, s0 );  
U = S(:,1);
V = S(:,2);
 
%   the linear stability analysis
p = (rho+kv)^2/rho;
u0 = (1-sqrt(1-4*p))/2;
v0 = sqrt(rho*(1-u0)/u0);
fu = -v0^2-rho;
fv = -2*u0*v0;
gu = v0^2;
gv= 2*u0*v0-rho-kv;

det = fu*gv-fv*gu  %this should be positive

us1 = (1+sqrt(1-4*p))/2;
vs1 = sqrt(rho*(1-us1)/us1);

figure(1)
 plot(u,v1,'--',u,v2,'--',U,V, 'linewidth',2)
 hold on
 plot(u0,v0,'*',us1,vs1,'*','linewidth',4)
axis([0 1 0 1])
xlabel('u','fontsize',20)
ylabel('v','fontsize',20)
  annotation('arrow',[ 0.65,0.57],[.45,0.53])
  
  stop
  
%Now plot the neutral stability (Turing) curve;

Dv = [0:.01:0.99]*gv;
Du = (det-Dv*fu)./(gv-Dv);

a = gv^2;
b = 2*fu*gv-4*det;
c = fu^2;

alph = (-b+sqrt(b^2-4*a*c))/(2*gv^2);

L = 20; % length of the domain
du = 1; % diffusion coefficient for u
dv =  0.025; % diffusion coefficient for v

figure(2)
 
plot(Dv,Du,Dv,alph*Dv,'--','linewidth',2)
hold on
axis([0 gv 0 4])
text(0.02,2,'unstable','fontsize',20)
text(0.06,1,'stable','fontsize',20)
xlabel('\omega^2D_v','fontsize',20)
ylabel('\omega^2D_u','fontsize',20)
set(gca,'Ytick',0:0.5:4)
for j = 1:15
    sc1=(j*pi/L)^2;
    plot(sc1*dv,sc1*du,'*','linewidth',4)
end
hold off
 
  
% now integrate the pde
 
N=150;  % number of spatial grid points
 
h = L/N;
    
 sc = [1;2*ones(N-2,1);1];
  X = h*(1:N)';
  %set initial data
  
  n0 = 8;  % set the wavelength
 U = u0-0.02*cos(n0*pi*X/L);
 V =  v0-0.02*cos(n0*pi*X/L);     
  %plot the solution at each time step

 figure(3)
  subplot(2,1,1)
  plot(X,V,'--','linewidth',2)
  xlabel('x','fontsize',20)
   ylabel('v(x,t)','fontsize',20)
  hold on
  subplot(2,1,2)
  plot(X, U,'--','linewidth',2)
   xlabel('x','fontsize',20)
    ylabel('u(x,t)','fontsize',20)
  axis([0 L  0 0.5])
 % axis([0 L 0 1])
 hold on 
 
tstep = 20; % time between plots
t_end = 2500; %total time to run simulation
 
%specify the output points
tspan = [0:tstep:t_end];

% integrate the system of ode's:
[T,S] = ode23(@pdeRHS,tspan,[U;V], odeset('maxstep',1));  

% plot the solution at each time step
 
for j = 1:length(T)
   figure(3)
  subplot(2,1,1)
  plot(X,S(j,N+1:2*N),'linewidth',2)
  xlabel('x','fontsize',20)
   ylabel('v(x,t)','fontsize',20)
  hold on
  subplot(2,1,2)
  plot(X, S(j,1:N),'linewidth',2)
   xlabel('x','fontsize',20)
    ylabel('u(x,t)','fontsize',20)
    pause(.1)
end
 
 figure(5) %plot the final solution
 plot(X, S(end,1:N),X,S(end,N+1:2*N),'linewidth',2)
 legend('u','v','fontsize',20)
 xlabel('x','fontsize',20)
 
end


%the right hand side for pde (MoL) simulation:
function s_prime=pdeRHS(t,s)
global   N sc rho kv L du dv 
 
h = L/N;
 
scu = du/h^2;
scv = dv/h^2;
U = s(1:N);
V = s(N+1:2*N);
F = U.*V.^2;
 Fu = scu*(-sc.*U+[0;U(1:end-1)]+[U(2:end);0])-F +rho*(1-U);
Fv = scv*(-sc.*V+[0;V(1:end-1)]+[V(2:end);0]) +F -(rho+kv)*V ;
s_prime = [Fu;Fv];
end

%the right hand side for ode simulation:
function s_prime=deRHS(t,s)
global rho kv
 
u = s(1);
v = s(2);
 
Fu = -u*v^2+rho*(1-u);
Fv = u*v^2-(rho+kv)*v;
 
 
s_prime = [Fu Fv]';

end
