% phase portrait for tiger bush ode
function tgpp
global a m N
set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 

 
%parameters
a = 6.5;
a=6;
m = 2;
 
w = [0.001:.01:a];
n1 = m./w;
w2 = [0.001:0.01:a];
n2 = sqrt(a./w2-1);

wp = (a+sqrt(a^2-4*m^2))/2;  
wm = (a-sqrt(a^2-4*m^2))/2;  

np = m/wp;
nm = m/wm;
% draw nullclines
figure(1)
plot(n1,w,'--',n2,w2,'--','linewidth',2)
hold on
plot( np,wp,'*',nm,wm,'*',...
    0,a,'*','linewidth',4)
axis([0 4 0 a])
 annotation('arrow',[ 0.55,0.65],[.45,0.37])
  annotation('arrow',[ 0.25,0.22],[.18,0.27])
 
xlabel('n','fontsize',20)
ylabel('w','fontsize',20)
%title('Tiger Bush Phase Plane','fontsize',20)
hold on

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

N = 2;  % number of trajectories
s0 = [1;1.4;0.01;0.01];
 
 
[T,S] = ode23(@deRHS,tspan, s0, odeset('maxstep',1));  

% plot trajectories
figure(1)
plot(S(:,1),S(:,3),'k',S(:,2),S(:,4),'k','linewidth',2)
hold off
 
 
%stability diagram
gm = [2:.01:10];
am = [0:.01:10];
ma = [2:.01:10];
atr = ma.^2./(sqrt(ma-1));  %critical stability curve

figure(2)
plot(am,am/2,atr,ma,'--', 'linewidth',2)
axis([0 10 0 5])
text(3,3,'I','fontsize',20)
%text(2,3,'no vegetation','fontsize',18)
text(6,1.5,'II','fontsize',20)
%text(4,1.5,'stable vegetation','fontsize',18)
text( 9,4.25,'III','fontsize',20)
%text(4,1.5,'unstable homogeneous vegetation','fontsize',18)
xlabel('a','fontsize',20)
ylabel('m','fontsize',20)
 
 
 

function fprime=deRHS(t,x)
global a m N
n = x(1:N);
w = x(N+1:2*N);

Fn = w.*n.^2-m*n;
Fw = a-w-w.*n.^2;

fprime = [Fn;Fw];
