% cell polarization ode system
clear all

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 

 
%parameters
del = 1;  % wlog
K = 1; % wlog
s = 0.06;

% plot the right hand side of the ode for several parameter values
 gam =  3;
Sset=[0.01,0.05,0.08];
Vt = 1.3;
for j = 1:length(Sset)
       St = Sset(j);
a = [0:.001:Vt-.01];
f = gam*a.^2./(K^2+a.^2);
g = (St+f).*(Vt-a)-a;
figure(1)
plot(a,g,'linewidth',2)
hold on
xlabel('u','fontsize',20)
ylabel('f(u,V_T-u)','fontsize',20)
end
plot([0,Vt],[0,0],'--')
legend('S=0.01','S=0.05','S=0.08','fontsize',16)
axis([0 0.8 -0.1 0.1])
 hold off
  

% %
% now plot steady state solution as a function of S
 Vtset = [  1.1;1.3;1.4];
 
 
for j = 1:length(Vtset)
       Vt = Vtset(j);
a = [0:.001:Vt-.01];
f = gam*a.^2./(K^2+a.^2);

  
S = -f+a./(Vt-a);
figure(2)
plot(S,a,'linewidth',2)
hold on
axis([0 0.2 0 1])
xlabel('S','fontsize',20)
ylabel('u','fontsize',20)
 end
legend('V_T=1.1','V_T=1.3','V_T=1.4','fontsize',16)
hold off
 

% now plot the critical curves separating different regions
s = 0.04;
 u = [10:95]/100/sqrt(3);

Vu = u.*(3-u.^2)./(1-3*u.^2);
gu = (s+1)*(1/3)*(3*u.^2-1).*(u.^2-3)./u.^2;
Vtu = 2*sqrt(gu+1)./gu;

figure(3)
%plot(gu,Vu, 'linewidth',2)
 plot( gu,Vtu,'linewidth',2)
xlabel('\gamma','fontsize',20)
ylabel('V_T','fontsize',20)
axis([0 4 0 4])
annotation('textarrow',[0.35 .41],[0.43 .53],'linewidth',2)
text(0.8,0.8,'I','fontsize',18)

text(1,1.4,'II','fontsize',18)
text(2.5,2.5,'III','fontsize',18)
 
% plot the Sp(u) double zero curve
gm=[0:0.1:4];
bt = - 8 *gm.^3 + 6* gm.^2 - 96*gm - 64;
al=27*gm.^2;
ct=27*gm.^2;

vsqp=max(0,( -bt+sqrt(max(0,bt.^2-4*al.*ct)))./(2*al));
vsqm=max(0,( -bt-sqrt(max(0,bt.^2-4*al.*ct)))./(2*al));
hold on
 
plot(gm,sqrt(vsqp),'linewidth',2)
%plot(gm,sqrt(vsqm),'--')
hold off

 


