% dd chemotaxis ss  phase plane
function chemotax
global K

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

 

ustar = 0.2;
K = log(ustar)/(ustar-1)
 u = [0:.01:2];
 f = exp(K*(u-1))-u;
 figure(1)
 plot(u,f)
 axis([0 2 -1 1])
 
tstep = 0.01;
t_end = 20;
tspan = [0:tstep:t_end];
N = 51;
dw=(1-ustar)/N;
w0 = ustar; 
L = zeros(N-1,1);
V = zeros(N-1,1);
for j = 1:N-1
   w0 = w0 + dw; 
s = [w0;0];

[T,S] = ode23s(@deRHS,tspan,s,odeset('Events',@events)); 
L(j) = T(end);
V(j) = S(end,1);
W(j) = w0;

 if (10*fix((j-1)/10)==j-1)
 figure(2)
 plot(S(:,1),S(:,2),S(:,1),-S(:,2),'linewidth',2)
 xlabel('v','fontsize',20)
 ylabel('dv/d\xi','fontsize',20)
 axis([0 1.5 -0.5 0.5])
 set(gca,'Xtick',0:0.25:1.5)
 set(gca,'Ytick',-0.5:0.25:0.5)
 hold on
 end

 if(j==18)
     figure(4)
 plot([T;2*T(end)-T([end:-1:1])], [S(:,1);S([end:-1:1],1)],...
     [T;2*T(end)-T([end:-1:1])],[exp(K*(S(:,1)-1));exp(K*(S([end:-1:1],1)-1))],'--','linewidth',2)
 xlabel('\xi','fontsize',20)
 ylabel('v, u/u^*','fontsize',20)
 legend('v','u/u^*','fontsize',20)
 set(gca,'Xtick',0:1:8)
 set(gca,'Ytick',0:0.25:2)
  
 end    
end

figure(3)
plot(2*L,V, 2*L,W,'linewidth',2)
xlabel('length l','fontsize',20)
 ylabel('v_{min},v_{max}','fontsize',20)
 axis([0 15 0 1.5])
set(gca,'Xtick',0:2.5:15)
 set(gca,'Ytick',0:0.25:1.5)
  
 
  figure(2)
  plot([0,2],[0,0],'--','linewidth',2)
plot([ustar,1],[0,0],'*','linewidth',2)
plot( [ustar,ustar],[-0.5,0.5],'--','linewidth',2)  
plot( [1,1],[-0.5,0.5],'--','linewidth',2) 
annotation('arrow',[.7,0.6],[0.17,0.17])
   annotation('arrow',[.6,0.7],[.88,0.88])
 hold off
function s_prime=deRHS(t,s)
global K R
u = s(1);
w = s(2);
 
Fu = w;
Fw = -exp(K*(u-1))+u;

s_prime = [Fu Fw]';


% event function
function [value,isterminal,direction]=events(t,y)

value = y(2);
isterminal = 1;
direction = 0;
