% this codes makes the first several of the plots in the Background chapter

function make_de_plots
global A

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


 
% plot the rhs of the de u' = f(u) 
x=[-.2:0.002:1.2];  % specify the x- coordinate values
N = length(x);
al = 0.25; %specify the parameter alpha
A = 10; % specify the scale factor
f = A*x.*(1-x).*(x-al); %the rhs function
figure(1)
plot(x,f, x,zeros(N,1),'--',[0,al,1],zeros(3,1),'*k','linewidth',2)

axis([-0.2 1.2 -1 1])
set(gca, 'XTick', -0.2 : 0.2 : 1.2);
xlabel('u','fontsize',20)
ylabel('du/dt','fontsize',20)
annotation('textarrow',[0.35 .25],[0.7 0.7],'linewidth',2)
annotation('textarrow',[0.55 .65],[0.7 0.7],'linewidth',2)

 
% the integral of 1/f(x):
F=(al*log(abs(1-x))-log(abs(x-al))+(1-al)*log(abs(x)))/(A*al*(al-1));
%plot F vs x
figure(2)
plot(x,F,x,zeros(N,1),'--','linewidth',2)

xlabel('u','fontsize',20)
ylabel('F(u)','fontsize',20)
set(gca, 'XTick', -0.2 : 0.2 : 1.2);
axis([-0.2 1.2 -3 2 ])
 
% reverse the axes:
figure(3)
 
plot(F,x,[-3,3],[al,al],'--',[-3,3],[0,0],'--',[-3,3],[1,1],'--','linewidth',2)

xlabel('t','fontsize',20)
ylabel('u(t)','fontsize',20)
axis([-3 2 -0.2 1.2])
 
 
% nullclines and flow directions for the de u'' + f(u) =0
figure(4)
plot(x,zeros(N,1),'--', [0,0],[-1,1],'--',[al,al],[-1,1],'--',[1,1],[-1,1],'--','linewidth',2)
hold on
plot([0,al,1],[0,0,0],'*','linewidth',4)
axis([-0.2 1.2 -1 1])
set(gca, 'XTick', -0.2 : 0.2 : 1.2);

annotation('textarrow',[0.215 0.265],[0.72 0.72],'linewidth',2)
annotation('textarrow',[0.36 0.41],[0.72 0.72],'linewidth',2)
annotation('textarrow',[0.77 0.83],[0.72 0.72],'linewidth',2)

annotation('textarrow',[0.265 0.215],[0.32 0.32],'linewidth',2)
annotation('textarrow',[0.41 0.36],[0.32 0.32],'linewidth',2)
annotation('textarrow',[0.83 0.77],[0.32 0.32],'linewidth',2)

annotation('textarrow',[0.6 .6],[0.56 0.49],'linewidth',2)
annotation('textarrow',[0.3 .3],[0.49 0.56],'linewidth',2)
annotation('textarrow',[0.18 .18],[0.56 0.49],'linewidth',2)
annotation('textarrow',[0.85 .85],[0.49 0.56],'linewidth',2)

annotation('textarrow',[0.28 .32],[0.6 0.63],'linewidth',2)
annotation('textarrow',[0.32 0.28],[0.39 0.42],'linewidth',2)
annotation('textarrow',[0.62 .58],[0.42 0.39],'linewidth',2)
annotation('textarrow',[0.58 .62],[0.63 0.6],'linewidth',2)
ylabel('v','fontsize',20)
xlabel('u','fontsize',20)
hold off
% make a cmyk file
 

%integral curves
F = (A/12)*x.^2.*(4*al*x-3*x.^2-6*al+4*x);
% create 4 curves

xj = [0,0.15,0.6,1]; % specify  initial points
for j = 1:4
    xs = xj(j);
    Fj = (A/12)*xs.^2.*(4*al*xs-3*xs.^2-6*al+4*xs);
    Fdf = Fj-F;
     
    ndx = find(Fdf>=0);
 ndxs = [ndx(2:end),0] ;
    xst = find(ndxs-ndx>1)
% plot intergral curves = trajectories    
figure(5)
if(isempty(xst))
    plot(x(ndx),sqrt(Fj-F(ndx)),x(ndx),-sqrt(Fj-F(ndx)),'linewidth',2)
else
    plot(x(ndx(1:xst)),sqrt(Fj-F(ndx(1:xst))),x(ndx(1:xst)),-sqrt(Fj-F(ndx(1:xst))),'linewidth',2)
plot(x(ndx( xst+1:end)),sqrt(Fj-F(ndx( xst+1:end))),x(ndx( xst+1:end)),-sqrt(Fj-F(ndx( xst+1:end))),'linewidth',2)
end
hold on
end
plot(x,zeros(N,1),'--', [0,0],[-1,1],'--',[al,al],[-1,1],'--',[1,1],[-1,1],'--','linewidth',2)
 annotation('textarrow',[0.53 .48],[0.36 0.33],'linewidth',2)
 annotation('textarrow',[0.48 .53],[0.72 0.68],'linewidth',2)
plot([0,al,1],[0,0,0],'*','linewidth',4)
axis([-0.2 1.2 -1 1])
set(gca, 'XTick', -0.2 : 0.2 : 1.2);
hold off
xlabel('u','fontsize',20)
ylabel('v','fontsize',20)
hold off
 


% now create some phase portraits
tstep = 0.01; %time step size
t_end = 5; % length of time interval
tspan = [0:tstep:t_end];

%saddle point
A = [0,1;1,0]  % the matrix rhs
eig(A)
s = [-1.13;1.1]; % initial data for integration
    
[T,S] = ode23s(@deRHS,tspan,s); 
 
figure(6)
plot(S(:,1),S(:,2),-S(:,1),-S(:,2),S(:,2),S(:,1),-S(:,2),-S(:,1),'linewidth',2)
 xlabel('u','fontsize',20)
ylabel('v','fontsize',20)
hold on
% add the stable and unstable manifolds
plot([0, 1],[0,1],[0,1],[0,-1],[0,-1],[0,1],[0,-1],[0,-1],'linewidth',2)
plot([-5,5],[0,0],'--',[0,0],[-5,5],'--',0,0,'*','linewidth',2) % nullclines
axis([-1 1 -1 1])
annotation('textarrow',[0.3 .3],[0.56 0.49],'linewidth',2)
annotation('textarrow',[0.75 .75],[0.49 0.56],'linewidth',2)
annotation('textarrow',[0.55 .49],[0.3 0.3],'linewidth',2)
annotation('textarrow',[0.49 .55],[0.75 0.75],'linewidth',2)
annotation('textarrow',[0.59 .65],[0.6 0.66],'linewidth',3)
annotation('textarrow',[0.65 0.59 ],[0.39 0.45],'linewidth',3)
 annotation('textarrow',[0.39 .45],[0.66 0.6],'linewidth',3)
 annotation('textarrow',[0.45 0.39 ],[0.455 0.395 ],'linewidth',3)

tstep = 0.01;
t_end = 35;
tspan = [0:tstep:t_end];
% unstable spiral
A = [0.1,1;-1,0.1]% the matrix rhs
eig(A)
s = [0.02;0.02];  %initial data
    
[T,S] = ode23s(@deRHS,tspan,s); 
 
figure(7)
plot(S(:,1),S(:,2),'linewidth',2)
 xlabel('u','fontsize',20)
ylabel('v','fontsize',20)
hold on
plot([-10,10],[1,-1],'--',[-1,1],[-10,10],'--',0,0,'*','linewidth',2) % nullclines
hold off
axis([-1 1 -1 1])
 
annotation('textarrow',[0.25 .25],[0.53 0.59],'linewidth',2)
annotation('textarrow',[0.75 .75],[0.54 0.47],'linewidth',2)
annotation('textarrow',[0.52 .46],[0.3 0.3],'linewidth',2)
annotation('textarrow',[0.51 .57],[0.7 0.7],'linewidth',2)
 

% stable spiral
A = [-0.1,1;-1,-0.1]
eig(A)

s = [-.4;-1];  %initial data
    
[T,S] = ode23s(@deRHS,tspan,s); 
 
figure(8)
plot(S(:,1),S(:,2),'linewidth',2)
 xlabel('u','fontsize',20)
ylabel('v','fontsize',20)
hold on
plot([-10,10],[-1,1],'--',[-1,1],[10,-10],'--',0,0,'*','linewidth',2) %nullclines
hold off
axis([-1 1 -1 1])
annotation('textarrow',[0.58 .52],[0.22 0.22],'linewidth',2)
annotation('textarrow',[0.47 .53],[0.78 0.78],'linewidth',2)
annotation('textarrow',[0.25 .25],[0.47 0.54],'linewidth',2)
annotation('textarrow',[0.73 .73],[0.58 0.51],'linewidth',2)
 
% stable node
A = [-1,0.2;0.3,-0.3]% the matrix rhs
eig(A)
idset = [0.4,-1,1,-1;-1,0,0,1];
 for j = 1:4
     s = idset(:,j); %specify initial data
[T,S] = ode23s(@deRHS,tspan,s); 
 
figure(9)
plot(S(:,1),S(:,2),'linewidth',2)
xlabel('u','fontsize',20)
ylabel('v','fontsize',20) 
hold on
 end

 % add nullclines and arrows
plot([-1,1],[-5,5],'--',[-10,10],[-10,10],'--',0,0,'*','linewidth',2)  % nullclines
hold off
axis([-1 1 -1 1])
annotation('textarrow',[0.46 .46],[0.19 0.27],'linewidth',2)
annotation('textarrow',[0.57 .57],[0.84 0.76],'linewidth',2)
annotation('textarrow',[0.27 .33],[0.3 0.3],'linewidth',2)
annotation('textarrow',[0.72 .65],[0.7 0.7],'linewidth',2)
 
function s_prime=deRHS(t,s)  % right hand side for ode system
global A
 
s_prime = A*s;