%% Appendix B. 100 more examples

%%
% <latex>
% \setcounter{page}{287}
% </latex>

ODEformats, fs = 9; lw = 0.7; ms = 6;

%%
% <latex>
% {\em 1. Van der Pol equation} [\kern .2pt Chap.~1\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(y) 0.05*diff(y,2)-(1-y^2)*diff(y)+y,[0 5]);
N.lbc = [1;0]; y = N\0; plot(y,LW,lw,CO,ivpnl)
axis square, hold off, set(gca,FS,6)
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,YT,-2:2)

%%
% <latex>
% {\em 2. Chirp} [\kern .2pt Chap.~2\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(t,y) diff(y)-cos(t^2)*y,[0 10]);
L.lbc = 1; y = L\0; plot(y,LW,lw,CO,ivp)
axis square, hold off, set(gca,FS,6)
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,YT,1:.5:3)

%%
% <latex>
% {\em 3. Growth with diminishing rate} [\kern .2pt Chap.~2\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(t,y) diff(y)-exp(-t/4)*y,[0 40]);
L.lbc = 1; y = L\0; plot(y,LW,lw,CO,ivp)
axis square, hold off, set(gca,FS,6)
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,YT,0:20:60)

%%
% <latex>
% {\em 4. Riccati equation} [\kern .2pt Chap.~3\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(t,y) diff(y)-y+t*y^2-t,[0 6]); N.lbc = 0;
y = N\0; plot(y,LW,1.2*lw,CO,ivpnl)
axis square, hold off, set(gca,FS,6)
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,YT,0:.5:1)

%%
% <latex>
% {\em 5. Bunching oscillations} (Bender \& Orszag p.~198) [\kern .2pt Chap.~3\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(t,y) diff(y)-cos(pi*t*y),[0 5]);
for a = 0:.2:2.8
  N.lbc = a; plot(N\0,LW,.8*lw,CO,ivpnl), hold on
end
axis square, hold off, set(gca,FS,6), ylim([0 3])
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,XT,0:5,YT,0:3)

%%
% <latex>
% {\em 6. Logistic equation} [\kern .2pt Chap.~3\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(y) diff(y)-(1-y)*y,[0 5]);
for a = .25:.25:2
  N.lbc = a; y = N\0; plot(y,LW,lw,CO,ivpnl), hold on
end
axis square, hold off, set(gca,FS,6), ylim([0 2])
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,XT,0:5,YT,0:.5:2)

%%
% <latex>
% {\em 7. Surprising complexity} [\kern .2pt Chap.~3\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(t,y) diff(y)-cos(t*y)*y-cos(t),[0 15]);
N.lbc = 0; y = N\0; plot(diff(y),CO,ivpnl,LW,lw)
axis square, hold off, set(gca,FS,6)
xlabel('$t$',FS,fs,IN,LT), ylabel('$y''$',FS,fs,IN,LT)
set(gca,YT,-1:2)

%%
% <latex>
% {\em 8. Equation with two steady states} (Davis, p.~$30$) [\kern .2pt Chap.~3\kern .3pt]
% \vskip -1.02em
% </latex>
d = [0 2.7]; N = chebop(@(t,y) diff(y)-t*y*(y-2),d);
t = chebfun('t',d);
for a = -2.5:.25:2
  N.lbc = a; y = N\0;
  plot(-flipud(t),flipud(y),LW,lw,CO,ivpnl), hold on
  arrowplot(t,y+1e-10,LW,lw,CO,ivpnl,MS,4)
end
axis square, hold off, set(gca,FS,6)
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 9. Equation with three steady states} [\kern .2pt Chap.~3\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(0,4);
N.op = @(t,y) diff(y)-y*abs(y)+t*y;
N.maxnorm = 4;
for y0 = -4:.1:4
  N.lbc = y0; y = N\0;
  plot(y,LW,.5*lw,CO,ivpnl), hold on, grid on
end
axis square, hold off, set(gca,FS,6), axis([0 4 -4 4])
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,XT,0:4,YT,-4:2:4)

%%
% <latex>
% {\em 10. Chase around the circle at 3/4 speed} [\kern .2pt Chap.~3\kern .3pt]
% \vskip -1.02em
% </latex>
d = [0 1.95*pi];
a = chebfun('exp(1i*t)',d);
N = chebop(@(t,z) diff(z)-.75*(a(t)-z)/abs(a(t)-z),d);
N.lbc = 0; arrowplot(a,LW,lw,CO,'r',YS,1.5), hold on
z = N\0; arrowplot(z,LW,lw,CO,ivpnl,YS,2)
axis square off, hold off, set(gca,FS,6)

%%
% <latex>
% {\em 11. Oscillation in a Lennard--Jones potential} [\kern .2pt Chap.~4\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(r) 0.1*diff(r,2)+(r/12)^(-1/12)-(r/12)^(-6),[0 12]);
N.lbc = [25;0]; r = N\0; plot(r,LW,1.2*lw,CO,ivp)
axis square, hold off, set(gca,FS,6), axis([0 12 0 30])
xlabel('$t$',FS,fs,IN,LT), ylabel('$r$',FS,fs,IN,LT)

%%
% <latex>
% {\em 12. Alternating growth and decay} [\kern .2pt Chap.~4\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(t,y) diff(y,2)+.06*sin(t/10)*diff(y)+y,[0 200]);
L.lbc = [1;0]; y = L\0; plot(y,LW,lw,CO,ivp)
axis square, hold off, set(gca,FS,6), ylim([-1.5 1.5])
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 13. Transient damping} [\kern .2pt Chap.~4\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(t,y) diff(y,2)+0.2*exp(-(t-50)^2/16)*diff(y)+y,[0 100]);
L.lbc = [1;0]; y = L\0; plot(y,LW,lw,CO,ivp)
axis square, hold off, set(gca,FS,6), ylim([-1.5 1.5])
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 14. Washing hanging on the line} [\kern .2pt Chap.~5\kern .3pt]
% \vskip -1.02em
% </latex>
d = [0 3];
L = chebop(@(y) diff(y,2),d,1,1);
x = chebfun('x',d); rhs = .1+0*x+(abs(x-2)<.1);
y = L\rhs; plot(y,LW,1.2*lw,CO,bvp)
axis square, axis([d .5 1]), set(gca,FS,6,XT,0:4,YT,0:.25:1)
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 15. Third-order problem} [\kern .2pt Chap.~5\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(y) 0.5*diff(y,3)+diff(y,2)+diff(y)+y,[-1 1]);
L.lbc = [0;0]; L.rbc = 1;
y = L\0; plot(y,LW,lw,CO,bvp)
axis square, set(gca,FS,6,XT,-1:.5:1,YT,0:.25:1)
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 16. Gaussian} [\kern .2pt Chap.~5\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(x,y) diff(y,2)+x*diff(y)+y,[-2 2],0,1);
y = L\0; plot(y,LW,lw,CO,bvp)
axis square, hold off, set(gca,FS,6,XT,-2:2,YT,0:4)
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 17. Beam with four interpolation points (= cubic spline)} [\kern .2pt Chap.~5\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(y) diff(y,4),[0 3]);
L.bc = @(x,y) [y(0); y(1); y(2)-1; y(3)-1];
y = L\0; plot(y,LW,lw,CO,bvp)
hold on, plot(0:3,[0 0 1 1],'.k',MS,12)
axis square, hold off, set(gca,FS,6,XT,0:3,YT,0:.5:1)
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 18. Troesch equation} [\kern .2pt Chap.~5\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(y) diff(y,2)-6*sinh(6*y),[0 1],0,1);
y = N\0; plot(y,LW,lw,CO,bvpnl)
axis square, hold off, set(gca,FS,6,XT,0:3,YT,0:.5:1)
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,XT,0:.25:1,YT,0:.25:1)

%%
% <latex>
% {\em 19. Harmonic oscillator for Schr\"odinger equation} [\kern .2pt Chap.~6\kern .3pt]
% \vskip -1.02em
% </latex>
x = chebfun('x',[-3 3]); V = x^2; h = 0.1;
L = chebop(@(y) -h^2*diff(y,2)+V*y,[-3 3]);
L.lbc = 0; L.rbc = 0; [W,D] = eigs(L,10);
for k = 1:10, plot(D(k,k)+0.06*W{k},LW,lw), hold on, end
plot(V,'k',LW,lw)
axis square, hold off, set(gca,FS,6), axis([-2 2 0 2.5])
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs-1,IN,LT)

%%
% <latex>
% {\em 20. V-shaped oscillator via {\tt quantumstates}} [\kern .2pt Chap.~6\kern .3pt]
% \vskip -1.02em
% </latex>
x = chebfun('x',[-3 3]); V = abs(x); quantumstates(V); title(' ')
axis square, hold off, set(gca,FS,6)
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 21. Lennard--Jones eigenstates} [\kern .2pt Chap.~6\kern .3pt]
% \vskip -1.02em
% </latex>
r = chebfun('r',[10 32]); V = (r/12)^(-12)-(r/12)^(-6); quantumstates(V); title(' ')
axis square, hold off, set(gca,FS,6), axis([10 30 -.28 .05])
xlabel('$r$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,YT,-.2:.1:0)

%%
% <latex>
% {\em 22. Orr-Sommerfeld eigenvalues} [\kern .2pt Chap.~6\kern .3pt]
% \vskip -1.02em
% </latex>
Re = 5772;           
B = chebop(-1,1); B.op = @(x,u) diff(u,2)-u; A = chebop(-1,1);
A.op = @(x,u) (diff(u,4)-2*diff(u,2)+u)/Re-1i*(2*u+(1-x^2)*(diff(u,2)-u));
A.lbc = [0;0]; A.rbc = [0;0]; lam = eigs(A,B,47);
plot(lam,'.r',MS,10), grid on, axis([-1 0 -1 0])
axis square, hold off, set(gca,FS,6)
xlabel('$\hbox{Re}(\lambda)$',FS,fs,IN,LT)
ylabel('$\hbox{Im}(\lambda)$',FS,fs,IN,LT)

%%
% <latex>
% {\em 23. Bessel equation} [\kern .2pt Chap.~7\kern .3pt]
% \vskip -1.02em
% </latex>
nu = 50;
L = chebop(@(x,y) x^2*diff(y,2)+x*diff(y)+(x^2-nu^2)*y,[0 150],0,1);
y = L\0; plot(y,LW,lw,CO,bvp)
axis square, hold off, set(gca,FS,6), ylim([-8 8])
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 24. Airy-like equation} (Lee \& Greengard) [\kern .2pt Chap.~7\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(x,y) 1e-4*diff(y,2)+(x^2-0.25)*y,[-1,1]);
L.lbc = 1; L.rbc = 2; y = L\0; plot(y,LW,.6*lw,CO,bvp)
axis square, hold off, set(gca,FS,6), ylim([-10 10])
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 25. Nearly ill-posed problem} [\kern .2pt Chap.~7\kern .3pt]
% \vskip -1.02em
% </latex>
x = chebfun('x',[-1 1]);
L = chebop(@(x,y) -.05*diff(y,2)+x*diff(y)-y,[-1 1]);
f = chebfun('exp(x)'); L.lbc = 0; L.rbc = 0; y = L\f;
plot(y,LW,lw,CO,bvp)
axis square, hold off, set(gca,FS,6), ylim([-400 400])
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,YT,-400:200:400)
%% 

%%
% <latex>
% {\em 26. Forcing at the resonant frequency} [\kern .2pt Chap.~8\kern .3pt]
% \vskip -1.02em
% </latex>
d = [0 80]; t = chebfun('t',d);
L = chebop(@(y) diff(y,2)+y,d); L.lbc = [0;0];
y = L\sin(t); plot(y,LW,lw,CO,ivp)
axis square, hold off, set(gca,FS,6), ylim([-40 40])
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,XT,0:20:80,YT,-40:20:40)

%%
% <latex>
% {\em 27. The same with damping} [\kern .2pt Chap.~8\kern .3pt]
% \vskip -1.02em
% </latex>
d = [0 80]; t = chebfun('t',d);
L = chebop(@(y) diff(y,2)+.08*diff(y)+y,d); L.lbc = [0;0];
y = L\sin(t); plot(y,LW,lw,CO,ivp)
axis square, hold off, set(gca,FS,6), ylim([-40 40])
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,XT,0:20:80,YT,-40:20:40)

%%
% <latex>
% {\em 28. Duffing equation (Davis p.~$399$)} [\kern .2pt Chap.~8\kern .3pt]
% \vskip -1.02em
% </latex>
d = [0 30]; t = chebfun('t',d);
N = chebop(@(y) diff(y,2)+y-.088*y^3,d); N.lbc = [-1;0];
y = N\sin(0.05*t); arrowplot(y,diff(y),LW,lw,CO,ivpnl,MS,5,YS,2.5)
axis square, hold off, set(gca,FS,6), axis([-1.5 5.5 -2 5])
xlabel('$y$',FS,fs,IN,LT), ylabel('$y''$',FS,fs,IN,LT)
%%

%%
% <latex>
% {\em 29. Forcing at the resonant frequency (phase plane)} [\kern .2pt Chap.~9\kern .3pt]
% \vskip -1.02em
% </latex>
d = [0 80]; t = chebfun('t',d);
L = chebop(@(y) diff(y,2)+y,d); L.lbc = [0;0];
y = L\sin(t); plot(y,diff(y),LW,lw,CO,ivp), axis(40*[-1 1 -1 1])
axis square, hold off, set(gca,FS,6)
xlabel('$y$',FS,fs,IN,LT), ylabel('$y''$',FS,fs,IN,LT)
set(gca,YT,-40:20:40)

%%
% <latex>
% {\em 30. The same with damping} [\kern .2pt Chap.~9\kern .3pt]
% \vskip -1.02em
% </latex>
d = [0 80]; t = chebfun('t',d);
L = chebop(@(y) diff(y,2)+.08*diff(y)+y,d); L.lbc = [0;0];
y = L\sin(t); plot(y,diff(y),LW,.8*lw,CO,ivp), axis(13*[-1 1 -1 1])
axis square, hold off, set(gca,FS,6)
xlabel('$y$',FS,fs,IN,LT), ylabel('$y''$',FS,fs,IN,LT)

%%
% <latex>
% {\em 31. Van der Pol equation limit cycle} [\kern .2pt Chap.~9\kern .3pt]
% \vskip -1.02em
% </latex>
T = 16;
N = chebop(@(y) diff(y,2)-4*(1-y^2)*diff(y)+y,[0 T]);
N.lbc = [1;1];
y = N\0; y = y{4,T}; plot(y,diff(y),LW,lw,CO,bvpnl)
hold on, plot(0,0,'.k',MS,ms)
axis square, hold off, set(gca,FS,6), axis([-7 7 -7 7])
xlabel('$y$',FS,fs,IN,LT), ylabel('$y''$',FS,fs,IN,LT)

%%
% <latex>
% {\em 32. Van der Pol equation in the stiff regime} [\kern .2pt Chap.~9\kern .3pt]
% \vskip -1.02em
% </latex>
T = 4;
N = chebop(@(y) .001*diff(y,2)-(1-y^2)*diff(y)+y,[0 T]);
N.lbc = [2;0];
pref = cheboppref; pref.ivpSolver = 'ode15s';
y = solveivp(N,0,pref); plot(y,LW,lw,CO,ivpnl,'jumpline','-')
axis square, hold off, set(gca,FS,6)
xlabel('$y$',FS,fs,IN,LT), ylabel('$y''$',FS,fs,IN,LT)

%%
% <latex>
% {\em 33. Rayleigh equation limit cycle} [\kern .2pt Chap.~9\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(y) diff(y,2)+6*(diff(y)^3-diff(y))+y,[0 20]);
N.maxnorm = 20;
N.lbc = [3;1];
y = N\0; y = y{4,18}; plot(y,diff(y),LW,lw,CO,bvpnl)
hold on, plot(0,0,'.k',MS,ms)
axis square, hold off, set(gca,FS,6), axis(3.5*[-1 1 -1 1])
xlabel('$y$',FS,fs,IN,LT), ylabel('$y''$',FS,fs,IN,LT)
set(gca,XT,-3:3,YT,-3:3)

%%
% <latex>
% {\em 34. Nonlinear pendulum heteroclinic paths} [\kern .2pt Chap.~9\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(y) diff(y,2)+sin(y),[0 21]);
N.lbc = [-3.1;0];
y = N\0; arrowplot(y,diff(y),LW,lw,CO,bvpnl,MS,4,'multi',20,YS,3)
hold on, plot(0,0,'.k',MS,ms)
axis square, hold off, set(gca,FS,6), axis(3.5*[-1 1 -1 1])
xlabel('$y$',FS,fs,IN,LT), ylabel('$y''$',FS,fs,IN,LT)

%%
% <latex>
% {\em 35. Unforced Duffing equation, phase portrait} [\kern .2pt Chap.~9\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(y) diff(y,2)+y-.088*y^3,[0 20]);
plot(0,0,'.k',MS,ms), hold on
for a = 0.28:.3:2.6
  N.lbc = [0; a]; y = N\0;
  arrowplot(y,diff(y),LW,lw,CO,ivpnl,MS,3,YS,3)
end
axis square, hold off, set(gca,FS,6,YT,-4:2:4), axis([-4 4 -4 4])
xlabel('$y$',FS,fs,IN,LT), ylabel('$y''$',FS,fs,IN,LT)

%%
% <latex>
% {\em 36. Cubic oscillator} [\kern .2pt Chap.~9\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(t,y) diff(y,2)+y^3,[0 20]);
for a = 0.1:.2:1.5
  N.lbc = [0;a];
  y = N\0; arrowplot(y,diff(y),LW,.8*lw,CO,ivpnl,MS,3,YS,2), hold on
end
hold on, plot(0,0,'.k',MS,ms)
axis square, hold off, set(gca,FS,6), axis([-2 2 -2 2])
xlabel('$y$',FS,fs,IN,LT), ylabel('$y''$',FS,fs,IN,LT)
set(gca,YT,-2:2)

%%
% <latex>
% {\em 37. Close to a saddle point} (Seydel p.~82) [\kern .2pt Chap.~9\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(y) diff(y,2)-y+y^2+0.45*diff(y)-0.5*y*diff(y),[0 14.5]);
N.lbc = [0;.3]; y = N\0; arrowplot(y,diff(y),LW,lw,MS,5,CO,ivpnl,YS,4)
hold on, plot(0,0,'.k',MS,ms)
axis square, hold off, set(gca,FS,6), axis([-1 2 -1 1])
xlabel('$y$',FS,fs,IN,LT), ylabel('$y''$',FS,fs,IN,LT)

%%
% <latex>
% {\em 38. Stable spiral point (scalar second-order)} [\kern .2pt Chap.~9\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(y) diff(y,2)+.1*diff(y,1)+y,[0 29]);
L.lbc = [0;1]; y = L\0; arrowplot(y,diff(y),LW,lw,CO,ivp,MS,5,YS,1.6)
ylim([-1 1.1]), hold on, plot(0,0,'.k',MS,ms)
axis square, hold off, set(gca,FS,6), axis(1.1*[-1 1 -1 1])
xlabel('$y$',FS,fs,IN,LT), ylabel('$y''$',FS,fs,IN,LT)

%%
% <latex>
% {\em 39. Walk on the sphere} [\kern .2pt Chap.~10\kern .3pt]
% \vskip -1.02em
% </latex>
d = [0 60];
L = chebop(d); L.lbc = [1;0;0];
t = chebfun('t',d);
f = 1; g = 0; h = 2*sin(t);
L.op = @(t,x,y,z) [diff(x)-f*y-g*z; diff(y)+f*x-h*z; diff(z)+g*x+h*y];
[x,y,z] = L\0; plot3(x,y,z,CO,ivp,LW,.6*lw)
hold on, plot3(x([0 end]),y([0 end]),z([0 end]),'.r',MS,10)
axis equal off, hold off, set(gca,FS,8)

%%
% <latex>
% {\em 40. Stable spiral point (first-order system)} [\kern .2pt Chap.~10\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(t,u,v) [diff(u)-v; diff(v)+u+.1*v],[0 29]);
L.lbc = [0;1]; [u,v] = L\0;
arrowplot(u,v,LW,lw,CO,ivp,MS,5,YS,2)
ylim([-1 1.1]), hold on, plot(0,0,'.k',MS,ms)
axis square, hold off, set(gca,FS,6), axis(1.1*[-1 1 -1 1])
xlabel('$u$',FS,fs,IN,LT), ylabel('$v$',FS,fs,IN,LT)

%%
% <latex>
% {\em 41. Brusselator chemical reaction} [\kern .2pt Chap.~10\kern .3pt]
% \vskip -1.02em
% </latex>
a = 1; b = 3;
N = chebop(@(t,u,v) [diff(u)-a+(b+1)*u-u^2*v; diff(v)-b*u+u^2*v],[0 15.3]);
N.lbc = [1;3.07];
[u,v] = N\0; arrowplot(u,v,LW,lw,CO,ivpnl,MS,5,YS,2.2)
hold on, plot([0 1],[0 3],'.k',MS,ms)
axis square, hold off, set(gca,FS,6), axis([0 4 .5 5])
xlabel('$u$',FS,fs,IN,LT), ylabel('$v$',FS,fs,IN,LT)
set(gca,XT,0:4,YT,0:5)

%%
% <latex>
% {\em 42. Full Brusselator} (Hairer--N\o rsett--Wanner) [\kern .2pt Chap.~10\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(0,60);
N.op = @(t,u) [diff(u{1})-1-u{1}^2*u{2}+(u{3}+1)*u{1};
               diff(u{2})-u{1}*u{3}+u{1}^2*u{2};
               diff(u{3})-1.2+u{3}*u{1}];
N.lbc = [1;2;2];
u = N\0; plot3(u{1},u{2},u{3},LW,lw)
axis square, hold off, set(gca,FS,6), axis([.5 2 .5 2 .5 2])
xlabel('$u_1$',FS,fs,IN,LT), ylabel('$u_2$',FS,fs,IN,LT), zlabel('$u_3$',FS,fs,IN,LT),
grid on, campos([-2.5 -9 8])
set(gca,XT,.5:.5:2,YT,.5:.5:2,'ztick',.5:.5:2)

%%
% <latex>
% {\em 43. Field--Noyes Oregonator -- a stiff equation} (Danby) [\kern .2pt Chap.~10\kern .3pt]
% \vskip -1.02em
% </latex>
s = 20; r = 0.1610; q = 8.375e-6;
N = chebop(@(t,u,v,w) ...
  [diff(u)-s*(v-u*v+u-q*u^2); ...
   diff(v)-(-v-u*v+w)/s; ...
   diff(w)-r*(u-w)],[0 250]);
N.lbc = [0;800;16];
pref = cheboppref; pref.ivpSolver = 'ode15s';
[u,v,w] = solveivp(N,0,pref); semilogy([u v w],LW,lw)
set(gca,'YMinorGrid','off')
axis square, hold off, set(gca,FS,6), axis([0 250 1e-3 1e6])
xlabel('$t$',FS,fs,IN,LT), ylabel('$u,v,w$',FS,fs,IN,LT)

%%
% <latex>
% {\em 44. Phase plane for nonautonomous 1st-order equation} [\kern .2pt Chap.~10\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(s,t,y) [diff(t)-1; diff(y)-exp(-t)+2*y],[0 3]);
quiver(L,[-1.5 2.5 -.5 1],'k','xpts',10,'ypts',10), hold on
L.lbc = [-1;0]; [t,y] = L\0;
arrowplot(t,y,LW,1.5*lw,CO,'r',YS,7)
axis square, hold off, set(gca,FS,6)
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 45. Homoclinic orbit} [\kern .2pt Chap.~10\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(t,u,v) [diff(u)-u.^3+2.*u.*v.^2; ...
                     diff(v)-2.*u.^2.*v+v.^3],[0 35]);
N.lbc = [0.16;0.01];
[u,v] = N\0; arrowplot(v,u,LW,lw,'multi',6,YS,2.5,MS,5), hold on, plot(0,0,'.k',MS,ms)
hold on, plot(0,0,'.k',MS,ms)
axis square, hold off, set(gca,FS,6), axis([-.1 1.1 -.1 1.1])
xlabel('$u$',FS,fs,IN,LT), ylabel('$v$',FS,fs,IN,LT)
set(gca,XT,0:.5:1,YT,0:.5:1)

%%
% <latex>
% {\em 46. Cubic damping} (Duan p.~6) [\kern .2pt Chap.~10\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(t,u,v) [diff(u)-v; diff(v)+u+v^3],[0 40]);
L.lbc = [0;1]; [u,v] = L\0;
arrowplot(u,v,LW,.8*lw,CO,ivpnl,MS,4,YS,2,'multi',4)
ylim([-1 1.1]), hold on, plot(0,0,'.k',MS,ms)
axis square, hold off, set(gca,FS,6), axis([-.8 1.05 -.75 1.1])
xlabel('$u$',FS,fs,IN,LT), ylabel('$v$',FS,fs,IN,LT)

%%
% <latex>
% {\em 47. Oscillatory system} (Hirsch--Smale--Devaney p.~229) [\kern .2pt Chap.~10\kern .3pt]
% \vskip -1.02em
% </latex>
a = 10; b = 1;
N = chebop(@(t,u,v) [diff(u)-a+u+4*u*v/(1+u^2); ...
                     diff(v)-b*u*(1-v/(1+u^2))],[0 11.8]);
N.lbc = [2.2;5]; [u,v] = N\0;
arrowplot(u,v,LW,lw,CO,ivpnl,MS,4,YS,1.8)
hold on, plot(2,5,'.k',MS,ms)
axis square, hold off, set(gca,FS,6), axis([0 4.5 3.5 8])
xlabel('$u$',FS,fs,IN,LT), ylabel('$v$',FS,fs,IN,LT)
set(gca,XT,0:4,YT,4:8)

%%
% <latex>
% {\em 48. Square limit cycle} (Hirsch--Smale--Devaney p.~215) [\kern .2pt Chap.~10\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(t,u,v) [diff(u)+sin(u)*(.1*cos(u)+cos(v)); ...
                     diff(v)+sin(v)*(.1*cos(v)-cos(u))],[0,60]);
N.lbc = [1.5;1.5]; [u,v] = N\0;
arrowplot(u,v,LW,lw,CO,ivpnl,MS,5,YS,2)
hold on, plot(pi/2,pi/2,'.k',MS,ms)
axis square, hold off, set(gca,FS,6), axis([-.5 3.5 -.5 3.5])
xlabel('$u$',FS,fs,IN,LT), ylabel('$v$',FS,fs,IN,LT)
set(gca,YT,0:3)

%%
% <latex>
% {\em 49. Four fixed points} (Davis p.~345) [\kern .2pt Chap.~10\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(t,u,v) [diff(u)+v-u^2-v^2; diff(v)-u+2*u*v],[0 5]);
N.maxnorm = [1.7 3];
for v0 = .5+.147*(-10:10)
  N.lbc = [0;v0]; [u,v] = N\0;
  plot(u,v,LW,.7*lw,CO,ivpnl), hold on, plot(-u,v,LW,.7*lw,CO,ivpnl)
end
hold on, plot([0 0 .5 -.5],[0 1 .5 .5],'.k',MS,ms)
axis square, hold off, set(gca,FS,6), axis([-1.7 1.7 -1.2 2.2]) 
set(gca,YT,-1:2)

%%
% <latex>
% {\em 50. Synchronizing three fireflies} (Strogatz, Danby) [\kern .2pt Chap.~10\kern .3pt]
% \vskip -1.02em
% </latex>
c = .2;
N = chebop(@(t,u,v,w) ...
      [diff(u)-0.97-c*sin(v-u)-c*sin(w-u);
       diff(v)-1.00-c*sin(u-v)-c*sin(w-v);
       diff(w)-1.03-c*sin(u-w)-c*sin(v-w)],[0 40]);
N.lbc = [0;2;4];
[u,v,w] = N\0; plot([sin(u) sin(v) sin(w)],LW,lw,CO,ivpnl)
axis square, hold off, set(gca,FS,6), axis([0 40 -1.5 1.5])
xlabel('$t$',FS,fs,IN,LT), ylabel('$u,v,w$',FS,fs,IN,LT)
set(gca,XT,0:20:80,YT,-1:1)

%%
% <latex>
% {\em 51. Hamiltonian system} (Jordan--Smith, p.~79) [\kern .2pt Chap.~10\kern .3pt]
% \vskip -1.02em
N = chebop(0,2);
N.op = @(t,u,v) [diff(u)-v*(13-u^2-v^2); diff(v)-12+u*(13-u^2-v^2)];
for a = -4.5:5.5
  N.lbc = [a;.2];
  [u,v] = N\0; plot(u,v,LW,.8*lw,CO,ivpnl), hold on
end
hold on, plot([-4 1 3],[0 0 0],'.k',MS,ms)
axis square, hold off, set(gca,FS,6), axis([-7.5 6.5 -7 7])
xlabel('$u$',FS,fs,IN,LT), ylabel('$v$',FS,fs,IN,LT)
set(gca,XT,-6:3:6,YT,-3:3:6)

%%
% <latex>
% {\em 52. A center and a saddle} [\kern .2pt Chap.~10\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(t,u,v) [diff(u)+u+v^2,diff(v)-v-u^2],[0 10]);
N.maxnorm = [3;3];
for a = -2:.15:0
  N.lbc = [2;a]; [u,v] = N\0; plot(u,v,LW,.8*lw,CO,ivpnl), hold on
  arrowplot(u{0,.5},v{0,.5},LW,.8*lw,CO,ivpnl,YS,2.5,MS,4)
end
for a = -.81:.15:-.3
  N.lbc = [-1;a]; [u,v] = N\0; plot(u,v,LW,.8*lw,'r')
  arrowplot(u{0,.5},v{0,.5},LW,.8*lw,CO,'r',YS,2.5,MS,4)
end
hold on, plot([-1 0],[-1 0],'.k',MS,ms)
axis square, hold off, set(gca,FS,6), axis([-3 2 -3 2])
xlabel('$u$',FS,fs,IN,LT), ylabel('$v$',FS,fs,IN,LT)

%%
% <latex>
% {\em 53. Leaky bucket} [\kern .2pt Chap.~11\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(y) diff(y)+y/sqrt(abs(y)),[0 4]);
N.lbc = 1; y = N\0; plot(y,LW,1.2*lw,CO,ivpnl)
axis square, hold off, set(gca,FS,6), axis([0 4 -.2 1.2])
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,XT,0:4,YT,0:.25:1)

%%
% <latex>
% {\em 54. Smooth random walk on the sphere} [\kern .2pt Chap.~12\kern .3pt]
% \vskip -1.02em
% </latex>
d = [0 2];
rng(6), u0 = randn(3,1); u0 = u0/norm(u0);
L = chebop(d); L.lbc = u0;
R = randnfun(0.05,3,d,'norm'); f = R(:,1); g = R(:,2); h = R(:,3);
L.op = @(t,x,y,z) [diff(x)-f*y-g*z; diff(y)+f*x-h*z; diff(z)+g*x+h*y];
[x,y,z] = L\0; plot3(x,y,z,CO,ivp,LW,.9*lw)
hold on, plot3(x([0 end]),y([0 end]),z([0 end]),'.r',MS,10)
axis equal off, hold off, set(gca,FS,8)

%%
% <latex>
% {\em 55. Less smooth random walk on the sphere} [\kern .2pt Chap.~12\kern .3pt]
% \vskip -1.02em
% </latex>
d = [0 2];
rng(6), u0 = randn(3,1); u0 = u0/norm(u0);
L = chebop(d); L.lbc = u0;
R = randnfun(0.005,3,d,'norm'); f = R(:,1); g = R(:,2); h = R(:,3);
L.op = @(t,x,y,z) [diff(x)-f*y-g*z; diff(y)+f*x-h*z; diff(z)+g*x+h*y];
[x,y,z] = L\0; plot3(x,y,z,CO,ivp,LW,.6*lw)
hold on, plot3(x([0 2]),y([0 2]),z([0 2]),'.r',MS,10)
axis equal off, hold off, set(gca,FS,6)

%%
% <latex>
% {\em 56. Nonuniqueness and extinction} [\kern .2pt Chap.~11\kern .3pt]
% \vskip -1.02em
% </latex>
S = chebfun('cos(pi*x/2)^2',[0 1]);
for h = 0:.125:1, plot([h*S 1-h*(1-S)],LW,.8*lw,CO,ivpnl), hold on, end
axis square, hold off, set(gca,FS,6), axis([0 1 0 1])
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,XT,0:.25:1,YT,0:.25:1)

%%
% <latex>
% {\em 57. Smooth Brownian motion to the circle} [\kern .2pt Chap.~12\kern .3pt]
% \vskip -1.02em
% </latex>
d = [0 2]; lam = 1e-3;
rng(0), f = randnfun(lam,d,'complex','norm');
L = chebop(@(y) diff(y),d); L.lbc = 0; L.maxnorm = 1;
plot(chebfun('exp(1i*pi*x)'),'k',LW,lw), hold on
y = L\f; plot(y,LW,.6*lw,CO,ivp)
axis square, hold off, set(gca,FS,6), axis off

%%
% <latex>
% {\em 58. BVP with smooth random coefficient} [\kern .2pt Chap.~12\kern .3pt]
% \vskip -1.02em
% </latex>
rng(0), d = [0 10]; f = randnfun(1,d);
L = chebop(@(y) .005*diff(y,2)+f*y,d); L.lbc = 0; L.rbc = 0;
y = L\1; plot(y,LW,lw,CO,bvp)
axis square, hold off, set(gca,FS,6)
set(gca,XT,0:2.5:10,YT,-20:20:40)

%%
% <latex>
% {\em 59. Complex decay with sidewise randomness} [\kern .2pt Chap.~12\kern .3pt]
% \vskip -1.02em
% </latex>
d = [0 40]; lam = 0.1;
rng(0), f = randnfun(lam,d,'norm');
L = chebop(@(y) diff(y)+(.3+1i*f)*y,d); L.lbc = 1; 
plot(chebfun('exp(1i*pi*x)'),'k',LW,.6*lw), hold on
y = L\0; plot(y,LW,lw,CO,ivp)
axis square, hold off, set(gca,FS,6), axis off

%%
% <latex>
% {\em 60. Conveyor belt} [\kern .2pt Chap.~12\kern .3pt]
% \vskip -1.02em
% </latex>
d = [0 1]; lam = 0.005;
rng(0), f = randnfun(lam,d,'norm'); g = randnfun(lam,d,'norm');
N = chebop(@(t,u,v) [diff(u)-60*exp(-40*v^2); diff(v)],d);
N.lbc = [0;0];
[u,v] = N\[f;g]; plot(u,v,LW,.6*lw,CO,ivpnl)
axis square, hold off, set(gca,FS,6,YT,-2:2), axis([0 6 -.5 2])
xlabel('$u$',FS,fs,IN,LT), ylabel('$v$',FS,fs,IN,LT)
set(gca,XT,0:2:6,YT,-.5:.5:2)

%%
% <latex>
% {\em 61. Equation with random $\pm 1$ forcing} [\kern .2pt Chap.~12\kern .3pt]
% \vskip -1.02em
% </latex>
rng(0), f = sign(randnfun(.5,[0,8],'norm'));
L = chebop(@(t,y) diff(y)-f*y,[0 8]);
L.lbc = 1; y = L\0; plot(y,CO,ivp,LW,.8*lw)
axis square, hold off, set(gca,FS,6)

%%
% <latex>
% {\em 62. Equation with intermittent kicks} [\kern .2pt Chap.~12\kern .3pt]
% \vskip -1.02em
% </latex>
d = [0 10]; t = chebfun('t',d);
rng(2), f = randnfun(.1,d,'norm')/(1.1-cos(2*pi*t));
L = chebop(@(t,y) diff(y)-f,d);
L.lbc = 0; y = L\0; plot(y,CO,ivp,LW,.8*lw)
axis square, hold off, set(gca,FS,6)
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT), ylim([-8 8])
set(gca,YT,-6:3:6)

%%
% <latex>
% {\em 63. Smooth geometric Brownian motion} [\kern .2pt Chap.~12\kern .3pt]
% \vskip -1.02em
% </latex>
d = [0 10]; lam = 0.05;
rng(2), f = randnfun(lam,d,'norm');
L = chebop(@(y) diff(y)-f*y,d); L.lbc = 1;
y = L\0; plot(y,LW,.7*lw,CO,ivp)
axis square, hold off, set(gca,FS,6),ylim([-1 3])
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,YT,-1:3), axis auto

%%
% <latex>
% {\em 64. Scalar chaotic equation} (Elwakil--Salama--Kennedy) [\kern .2pt Chap.~13\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(t,y) 1.5*diff(y,3)+diff(y,2)+diff(y)+y-tanh(8*y),[0 200]);
N.lbc = [2;0;0];
y = N\0; plot(y,LW,lw,CO,ivpnl)
axis square, hold off, set(gca,FS,6), ylim([-2.5,2.5])
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,YT,-2:2)

%%
% <latex>
% {\em 65. Forced damped nonlinear pendulum} (Danby) [\kern .2pt Chap.~13\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(t,y) diff(y,2)+.1*diff(y)+sin(y)-sin(t),[0 200]);
N.lbc = [0; 0.74];
y = N\0; plot(y,LW,.8*lw,CO,ivpnl)
axis square, hold off, set(gca,FS,6), ylim([-10,20])
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,YT,-10:10:20)

%%
% <latex>
% {\em 66. R\"ossler equation period doubling} [\kern .2pt Chap.~13\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(0,250); N.lbc = [2;0;0];
N.op = @(t,u,v,w) [diff(u)+v+w; diff(v)-u-.2*v; diff(w)-.2-w*(u-3.5)];
[u,v,w] = N\0; plot(u{200,250},LW,lw,CO,ivpnl)
axis square, hold off, set(gca,FS,6), axis([200 250 -10 10])
xlabel('$t$',FS,fs,IN,LT), ylabel('$u$',FS,fs,IN,LT)
set(gca,XT,200:25:250,YT,-8:4:8)

%%
% <latex>
% {\em 67. 3-body problem} [\kern .2pt Chap.~13\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(0,4);
N.op = @(t,u,v,w) [ ...
    diff(u,2)+(u-v)/abs(u-v)^3+(u-w)/abs(u-w)^3; ...
    diff(v,2)+(v-u)/abs(v-u)^3+(v-w)/abs(v-w)^3; ...
    diff(w,2)+(w-u)/abs(w-u)^3+(w-v)/abs(w-v)^3];
N.lbc = @(u,v,w) [u-1; v+1; w-(0.1+1.5i); diff(u); diff(v); diff(w)];
[u,v,w] = N\0; plot([u;v;w],LW,lw)
axis square, hold off, set(gca,FS,6), axis([-1.7 3.2 -1.7 3.2])
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,XT,-1:3,YT,-1:3)

%%
% <latex>
% {\em 68. 3-body problem, slightly different data} [\kern .2pt Chap.~13\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(0,6);
N.op = @(t,u,v,w) [ ...
    diff(u,2)+(u-v)/abs(u-v)^3+(u-w)/abs(u-w)^3; ...
    diff(v,2)+(v-u)/abs(v-u)^3+(v-w)/abs(v-w)^3; ...
    diff(w,2)+(w-u)/abs(w-u)^3+(w-v)/abs(w-v)^3];
N.lbc = @(u,v,w) [u-1; v+1; w-(0.1+1i); diff(u); diff(v); diff(w)];
[u,v,w] = N\0; plot([u;v;w],LW,lw)
axis square, hold off, set(gca,FS,6), axis([-1.2 1.2 -0.8 1.6])
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,XT,-1:.5:1,YT,-1:.5:1)

%%
% <latex>
% {\em 69. Double pendulum} (Danby) [\kern .2pt Chap.~13\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(0,150);
c = @(u) cos(u); s = @(u) sin(u); cs = @(u) cos(u).*sin(u); d = @(u) diff(u);
N.op = @(t,u,v) ...
   [diff(u,2)+(d(u)^2*cs(u-v)+d(v)^2*s(u-v)+2*s(u)-c(u-v)*s(v))/(2-c(u-v)^2); ...
    diff(v,2)+(-2*d(u)^2*s(u-v)-d(v)^2*cs(u-v)-2*c(u-v)*s(u)+2*s(v))/(2-c(u-v)^2)];
N.lbc = @(u,v) [u-3.14159; v; diff(u); diff(v)]; [u,v] = N\0;
plot([u v],LW,lw)
axis square, hold off, set(gca,FS,6), axis([0 150 -17 17])
xlabel('$t$',FS,fs,IN,LT), ylabel('$u,v$',FS,fs,IN,LT)

%%
% <latex>
% {\em 70. H\'enon--Heiles equations} [\kern .2pt Chap.~13\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(t,u,v) [diff(u,2)+u+2*u*v; diff(v,2)+v+u^2-v^2],[0 100]);
N.lbc = @(u,v) [u-.32;v-.35;diff(u);diff(v)];
[u,v] = N\0; plot(u,v,LW,.6*lw,CO,ivpnl)
axis square, hold off, set(gca,FS,6), axis([-.68 .52 -.6 .6])
xlabel('$u$',FS,fs,IN,LT), ylabel('$v$',FS,fs,IN,LT)
set(gca,XT,-.6:.2:.6,YT,-.6:.2:.6)

%%
% <latex>
% {\em 71. Spiral sink} [\kern .2pt Chap.~15\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(t,u,v) [diff(u)+u/2-v; diff(v)+u-v/3],[0 30]);
L.lbc = [0;1];
[u,v] = L\0; arrowplot(u,v,LW,lw,MS,4,CO,ivp,YS,4)
hold on, plot(0,0,'.k',MS,8)
axis square, hold off, set(gca,FS,6), axis(1.3*[-1 1 -1 1])
xlabel('$u$',FS,fs,IN,LT), ylabel('$v$',FS,fs,IN,LT)
set(gca,XT,-1:.5:1,YT,-1:.5:1)
%%
% <latex>
% {\em 72. Spiral source} [\kern .2pt Chap.~15\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(t,u,v) [diff(u)-.4*u+v; diff(v)-u+.2*v],[0 30]);
L.lbc = [0;1];
[u,v] = L\0; arrowplot(u,v,LW,lw,MS,4,CO,ivp,YS,2.5)
hold on, plot(0,0,'.k',MS,8)
axis square, hold off, set(gca,FS,6), axis(23*[-1 1 -1 1])
xlabel('$u$',FS,fs,IN,LT), ylabel('$v$',FS,fs,IN,LT)
set(gca,XT,-20:10:20,YT,-20:10:20)

%%
% <latex>
% {\em 73. Linear transient growth} [\kern .2pt Chap.~15\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(t,u,v) [diff(u)+u-20*v; diff(v)+2*v],[0 3]);
L.lbc = [0;1];
[u,v] = L\0; arrowplot(u,v,LW,lw,MS,4,CO,ivp)
hold on, plot(0,0,'.k',MS,8)
axis square, hold off, set(gca,FS,6), axis([-1 6 -3.5 3.5])
xlabel('$u$',FS,fs,IN,LT), ylabel('$v$',FS,fs,IN,LT)
set(gca,XT,-1:6,YT,-3:3)

%%
% <latex>
% {\em 74. From slightly unstable to slightly stable} [\kern .2pt Chap.~15\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(t,y) diff(y,2)+.05*tanh(y)*diff(y)-sin(y),[0 70]);
N.lbc = [-pi;1]; y = N\0; arrowplot(y,diff(y),CO,ivpnl,YS,3,MS,4,LW,.8*lw)
hold on, plot(pi*(-1:1),0*(-1:1),'.k',MS,8)
axis square, hold off, set(gca,FS,6), axis([-6 6 -6 6])
xlabel('$y$',FS,fs,IN,LT), ylabel('$y''$',FS,fs,IN,LT)
%set(gca,XT,-1:6,YT,-3:3)

%%
% <latex>
% {\em 75. Euler buckling} [\kern .2pt Chap.~16\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(x,y) diff(y,2)+3*sin(y),[-1 1],0,0);
for k = 1:-1:-1
  N.init = chebfun(@(x) k*cos(pi*x/2)); plot(N\0,LW,1.2*lw), hold on
end
axis square, hold off, set(gca,FS,6), %ylim([-3 3])
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 76. Carrier equation} [\kern .2pt Chap.~16\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(x,y) 0.001*diff(y,2)+2*(1-x^2)*y+y^2,[-1 1]);
N.lbc = 0; N.rbc = 0; x = chebfun('x');
y = N\1; plot(y,LW,lw,CO,bvpnl)
axis square, hold off, set(gca,FS,6), ylim([-3 3])
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 77. Orbit around two fixed suns} [\kern .2pt Chap.~16\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(0,12);
N.op = @(t,x,y) ...
  [diff(x,2)+x/(x^2+y^2)^1.5+x/(x^2+(y-1)^2)^1.5; ...
   diff(y,2)+y/(x^2+y^2)^1.5+(y-1)/(x^2+(y-1)^2)^1.5];
N.lbc = @(x,y) [x+2; y-1; diff(x)-0.8; diff(y)-0.454]; hold off
[x,y] = N\0; subplot(1,2,1)
plot([0 0],[0 1],'.',CO,orange,MS,20), hold on, grid on
arrowplot(x,y,LW,lw,CO,ivpnl)
axis square, hold off, set(gca,FS,6), axis([-3.2 0.8 -0.8 3.2])
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,XT,-3:0,YT,0:3)

%%
% <latex>
% {\em 78. Hopf bifurcation linear} [\kern .2pt Chap.~17\kern .3pt]
% \vskip -1.02em
% </latex>
T = 100; dom = [-T,T];
L = chebop(@(t,y) diff(y,2)-.15*(t/T)*diff(y)+y,dom);
L.lbc = [0;0]; rng(0), f = 0.1*randnfun(dom,'norm');
y = L\f; plot(y,LW,lw,CO,ivp)
axis square, hold off, set(gca,FS,6), ylim([-40 40])
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,YT,-40:20:40)

%%
% <latex>
% {\em 79. Hopf bifurcation nonlinear system} (Seydel p.~76) [\kern .2pt Chap.~17\kern .3pt]
% \vskip -1.02em
% </latex>
T = 100; dom = [-T/5,T];
N = chebop(@(t,u,v) [diff(u)+v-u*(t/T-u^2-v^2); ...
                     diff(v)-u-v*(t/T-u^2-v^2)],dom);
N.lbc = [0;0]; rng(0), f = 0.01*randnfun(1,dom,'norm');
[u,v] = N\[f;0]; t = chebfun('t',dom);
plot3(t,u,v,LW,lw,CO,ivp), xlim([-20 100]), view(-12,44)
axis square, hold off, set(gca,FS,6)
xlabel('$t$',FS,fs,IN,LT), ylabel('$u$',FS,fs,IN,LT)
zlabel('$v$',FS,fs,IN,LT)
set(gca,XT,0:40:80,YT,-1:1,'ztick',-1:1)

%%
% <latex>
% {\em 80. Hopf bifurcation van der Pol eq.} [\kern .2pt Chap.~17\kern .3pt]
% \vskip -1.02em
% </latex>
T = 60; dom = [-T,T];
N = chebop(@(t,y) diff(y,2)-2*(t/T)*(1-y^2)*diff(y)+y,dom);
N.lbc = [0;0]; rng(0), f = 0.05*randnfun(dom,'norm');
y = N\f; plot(y,LW,lw,CO,ivpnl)
axis square, hold off, set(gca,FS,6)
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 81. Continuous stirred tank reactor} (Seydel p.~157) [\kern .2pt Chap.~17\kern .3pt]
% \vskip -1.02em
% </latex>
T = 50;
N = chebop(@(t,u,v) [diff(u)+u-(0.1+.2*t/T)*(1-u)*exp(v); ...
                     diff(v)+v-16.2*(0.1+.2*t/T)*(1-u)*exp(v)+3*v],[0 T]);
N.lbc = [.2;.6];
[u,v] = N\0; plot(u,LW,lw,CO,ivpnl), ylim([0.1 1.1])
axis square, hold off, set(gca,FS,6)
xlabel('$t$',FS,fs,IN,LT), ylabel('$u$',FS,fs,IN,LT)
set(gca,XT,0:10:50,YT,.2:.2:1)

%%
% <latex>
% {\em 82. Periodic variant of Fig.\ $7.5$} [\kern .2pt Chap.~19\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(x,y) 0.005*diff(y,2)-sin(x)*y,[-2*pi 2*pi]);
L.bc = 'periodic'; y = L\1; plot(y,LW,lw,CO,bvp), ylim([-40 40])
axis square, hold off, set(gca,FS,6)
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 83. Periodic ODE in polar coordinates} [\kern .2pt Chap.~19\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(theta,r) diff(r,2)+(r+.35)/(1.1+sin(6*theta)),[0 2*pi]);
theta = chebfun('t',[0 2*pi]);
L.bc = 'periodic'; r = L\1; plot(exp(1i*theta)*r,CO,bvp,LW,1.2*lw)
axis square, hold off, axis(.15*[-1 1 -1 1])
set(gca,FS,6,XT,-.1:.1:.1,YT,-.1:.1:.1)
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 84. Stable Hill equation} [\kern .2pt Chap.~19\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(t,y) diff(y,2)+exp(2.40*sin(t))*y,[0 100]);
L.lbc = [1;0]; y = L\0; plot(y,LW,.8*lw,CO,ivp)
axis square, hold off, set(gca,FS,6)
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 85. Unstable Hill equation} [\kern .2pt Chap.~19\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(t,y) diff(y,2)+exp(2.45*sin(t))*y,[0 100]);
L.lbc = [1;0]; y = L\0; plot(y,LW,.8*lw,CO,ivp)
axis square, hold off, set(gca,FS,6)
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 86. Periodic problem with discontinuous coefficient} [\kern .2pt Chap.~19\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(x,y) diff(y)+(abs(x-3)<1)*y,[0 2*pi]);
L.bc = 'periodic'; y = L\1; plot(y,CO,bvp,LW,1.1*lw)
axis square, hold off, set(gca,FS,6,YT,0:2:8), ylim([0 8])
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 87. Periodic Carrier equation} [\kern .2pt Chap.~19\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(x,y) diff(y,2)+2*(2-cos(x))*y+y^2,[-pi pi],'periodic');
N.init = chebfun('-cos(2*x)',[-pi pi],'trig');
y = N\1; plot(y,LW,lw,CO,bvpnl)
axis square, hold off, set(gca,FS,6), ylim(2.2*[-1 1])
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,YT,-2:2)

%%
% <latex>
% {\em 88. Periodic Carrier equation, 2nd soln.} [\kern .2pt Chap.~19\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(x,y) diff(y,2)+2*(2-cos(x))*y+y^2,[-pi pi],'periodic');
N.init = chebfun('1.5*cos(2*x)',[-pi pi],'trig');
y = N\1; plot(y,LW,lw,CO,bvpnl)
axis square, hold off, set(gca,FS,6), ylim(2.2*[-1 1])
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,YT,-2:2)

%%
% <latex>
% {\em 89. Periodic Carrier equation, 3rd soln.} [\kern .2pt Chap.~19\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(x,y) diff(y,2)+2*(2-cos(x))*y+y^2,[-pi pi],'periodic');
N.init = chebfun('cos(x)^2',[-pi pi],'trig');
y = N\1; plot(y,LW,lw,CO,bvpnl)
axis square, hold off, set(gca,FS,6), ylim(2.2*[-1 1])
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,YT,-2:2)

%%
% <latex>
% {\em 90. Boundary layer at left} (Hemker) [\kern .2pt Chap.~20\kern .3pt]
% \vskip -1.02em
% </latex>
ep = 0.05;
L = chebop(@(x,y) ep*diff(y,2)+diff(y)-(1+ep)*y,[-1 1]);
L.lbc = 1+exp(-2); L.rbc = 1+exp(-2*(1+ep)/ep); y = L\0; plot(y,LW,lw,CO,bvp)
axis square, hold off, set(gca,FS,6), %ylim([-2 2])
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,XT,-1:.5:1)

%%
% <latex>
% {\em 91. Linear shock} (Lee \& Greengard) [\kern .2pt Chap.~20\kern .3pt]
% \vskip -1.02em
% </latex>
ep = 1e-4;
L = chebop(@(x,y) ep*diff(y,2)+2*x*diff(y),[-1 0 1],-1,1);
y = L\0; plot(y,LW,lw,CO,bvp)
axis square, hold off, set(gca,FS,6), ylim([-2 2])
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,YT,-30:15:30)

%%
% <latex>
% {\em 92. Linear corner layer} (Hemker) [\kern .2pt Chap.~20\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(x,y) 0.001*diff(y,2)+x*diff(y)-y,[-1 1],.8,1.2);
y = L\0; plot(y,LW,lw,CO,bvpnl)
axis square, hold off, set(gca,FS,6), ylim([-1,2])
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 93. Nonlinear corner layer} (Hemker) [\kern .2pt Chap.~20\kern .3pt]
% \vskip -1.02em
% </latex>
N = chebop(@(y) 0.05*diff(y,2)+diff(y)^2-1,[-1 1],.8,1.2);
y = N\1; plot(y,LW,lw,CO,bvpnl)
axis square, hold off, set(gca,FS,6), ylim([-1,2])
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)

%%
% <latex>
% {\em 94. Cusp} (Hemker, Lee \& Greengard) [\kern .2pt Chap.~20\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(x,y) 1e-6*diff(y,2)+x*diff(y)-.5*y,[-1 0 1],1,2);
y = L\0; plot(y,LW,lw,CO,bvp)
axis square, hold off, set(gca,FS,6), ylim([-0.5 2.5])
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
%% 
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% {\em 95. Ill-conditioned equation} (Lee \& Greengard) [\kern .2pt Chap.~20\kern .3pt]
% \vskip -1.02em
% </latex>
L = chebop(@(x,y) 0.02*diff(y,2)-x*diff(y)+y,[-1 1],1,2);
y = L\0; plot(y,LW,lw,CO,bvp)
axis square, hold off, set(gca,FS,6), ylim([-1 2])
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
%% 
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% {\em 96. Shock layer} (Hemker) [\kern .2pt Chap.~20\kern .3pt]
% \vskip -1.02em
% </latex>
ep = 1e-6;
L = chebop(@(x,y) ep*diff(y,2)+x*diff(y),[-1 0 1],-2,0);
f = chebfun(@(x) -ep*pi^2*cos(pi*x)-pi*x*sin(pi*x));
y = L\f; plot(y,LW,lw,CO,bvp)
axis square, hold off, set(gca,FS,6), ylim([-2.5 2.5])
xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,XT,-1:.5:1,YT,-2:2)
%% 
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% {\em 97. Blowup equation with complex perturbation} (Herzog \& Mattingley 2015) [\kern .2pt Chap.~21\kern .3pt]
% \vskip -1.02em
% </latex>
d = [0 1.03]; N = chebop(@(y) diff(y)-y^2,d);
N.lbc = 1; rng(0); f = 0.01*randnfun(d,'norm','complex');
y = N\f; arrowplot(y,LW,lw,YS,1.2,CO,ivp)
axis square, hold off, set(gca,FS,6), axis([-70 70 -20 120])
xlabel('$\hbox{Re}(y)$',FS,fs,IN,LT), ylabel('$\hbox{Im}(y)$',FS,fs,IN,LT)
set(gca,XT,-80:40:80,YT,0:40:120)

%%
% <latex>
% {\em 98. KdV train of solitons} [\kern .2pt Chap.~22\kern .3pt]
% \vskip -1.02em
% </latex>
d = [-5,5]; x = chebfun('x',d);
S = spinop(d,[0 2]);
S.lin = @(u) -0.01*diff(u,3); S.nonlin = @(u) -0.03*diff(u.^2);
c = 1; soliton = 50*c*sech(sqrt(100*c)*x/2)^2;
S.init = .25*soliton;
u = spin(S,160,.005,'plot','off'); plot([S.init u],LW,lw)
axis square, hold off, set(gca,FS,6), ylim([-5 15])
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,XT,-5:2.5:5',YT,-5:5:15)

%%
% <latex>
% {\em 99. Kuramoto--Sivashinsky chaotic state} [\kern .2pt Chap.~22\kern .3pt]
% \vskip -1.02em
% </latex>
d = [-50,50]; x = chebfun('x',d);
S = spinop(d,[0 100]);
S.lin = @(u) -diff(u,2)-diff(u,4); S.nonlin = @(u) -0.5*diff(u.^2);
S.init = exp(-(x/5)^2);
u = spin(S,200,.1,'plot','off'); plot(u,LW,lw,'k')
axis square, hold off, set(gca,FS,6), ylim([-4 4])
xlabel('$t$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
set(gca,XT,-50:25:50,YT,-4:2:4)

%%
% <latex>
% {\em 100. Ginzburg--Landau PDE} [\kern .2pt Chap.~22\kern .3pt]
% \vskip -1.02em
% </latex>

%%
%   u = spin2('gl','plot','off');
%   plot(real(u)), view(0,90)
%   colormap(parula)
%   axis square, hold off, set(gca,FS,6)
%   xlabel('$x$',FS,fs,IN,LT), ylabel('$y$',FS,fs,IN,LT)
%   set(gca,XT,0:50:100,YT,0:50:100)
