% fisher equation steady state  phase plane
function fisher_ss
global d sm sp

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

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

N = 500; % number of trajectories to find

 %there are three cases
for cse = 1:3
switch cse

    case 1%Case 1: Fisher w/ Dirichlet
d = 0;
c = 0;
dw=sqrt(1.00301/3)/N;
sp = 0;
sm = 0;

    case 2
        %Case 2: Fisher w/ Robin condition
c = 0;
d = 0.2; % delta = 1/d
dw=sqrt(0.967/3)/N;
sp = d;
sm = d;

    case 3
% Case 3: 
d = 0.4; %delta^2 = 1/d;
c = 1.0;
sp = 2*d/(-c*d+sqrt(c^2*d^2+4));
sm = -2*d/(-c*d-sqrt(c^2*d^2+4));
dw=sqrt(0.837/3)/N;
end
tstep = 0.01;
t_end = 20;
tspan = [0:tstep:t_end];

w0 = 0; 
L = zeros(N,1);
V = zeros(N,1);
for j = 1:N
   w0 = w0 + dw;
   v0 = sp*w0;
s = [v0;w0]; % solve the differential equations with these initial data

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

if(j==100*fix(j/100)) % plot only certain of the trajectories to plot
 figure(3*cse)
 %plot(S(:,1),S(:,2),S(:,1),-S(:,2),'linewidth',2)
 plot(S(:,1),S(:,2),[0,sp],[0,1],'--',...
     [0,sm],[0,-1],'--','linewidth',2)
 xlabel('v','fontsize',20)
 ylabel('w','fontsize',20)
 axis([0 1 -0.6 0.6])
 annotation('textarrow',[0.45,0.55],[0.78 0.73],'linewidth',2)
 annotation('textarrow',[0.55,0.45],[0.32 0.27],'linewidth',2)
 hold on
end

end
plot([0,1],[0,0],'--','linewidth',2)
plot([0,1],[0,0],'*','linewidth',2)

if(cse==1)
figure(4)
plot(L,V,'linewidth',2)
xlabel('Y','fontsize',20)
 ylabel('V_{max}','fontsize',20)
 axis([0 10 0 1])
 hold on
end
if(cse==2)
    figure(4)
plot(L,V,'--','linewidth',2)
xlabel('Y','fontsize',20)
 ylabel('V_{max}','fontsize',20)
 axis([0 10 0 1])
 hold off
end
 
if(cse==3)
 figure(10)
 cc = [0:.01:1.99];
 Y = 2*atan((4+d*cc.^2)./(4*d-d*cc.^2))./sqrt(4-cc.^2)/(atan(1/d)); % the critical curve
 plot(cc,Y,'linewidth',2)
 xlabel('C','fontsize',20)
 ylabel('Y/Y_0(\delta)','fontsize',20)
 text(0.6,2,'survival','fontsize',20)
 text(1.4,0.5,'extinction','fontsize',20)
 axis([0 2 0 4])
end
end

if(cmykflg==1) 
    figure(3)
print('../../figs_c/chapt_6/fisher_ss_pp','-deps','-cmyk')
figure(4)
print('../../figs_c/chapt_6/v_vs_l','-deps','-cmyk')
figure(6)
print('../../figs_c/chapt_6/fshr_ss_pp_dne0','-deps','-cmyk')
figure(9)
print('../../figs_c/chapt_6/fshr_ss_pp_cdne0','-deps','-cmyk')
figure(10)
print('../../figs_c/chapt_6/mvng_niche_crit','-deps','-cmyk')

end

function s_prime=deRHS(t,s)
v = s(1);
w = s(2);
 
Fv = w;
Fw = -v*(1-v);

s_prime = [Fv Fw]';


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

value = sm*y(2)+ y(1);  % this is the equation for the ending condition
isterminal = 1;
direction = 0;
