%% 10. Systems of equations
% <latex>
% \setcounter{page}{114}
% </latex>

%%
% <latex>
% Up to this point, our discussion has concerned scalar ODE\kern .5pt s.\ \ In this
% chapter we begin to talk about systems of equations, which
% involve more than one dependent variable, or
% equivalently, a dependent variable that is a vector rather than
% a scalar.  
% </latex>

%%
% <latex>
% All our vector problems will be of first order.  A general first-order
% ODE can be written
% $$ {\bf y}'(t) = {\bf f}(t,{\bf y}) \eqno (10.1) $$
% for a linear or nonlinear function $\bf f$, where for each
% $t$, ${\bf y}(t)$ 
% represents a vector of dimension $n$,
% $$ {\bf y}(t) = \pmatrix{y_1^{}(t)\cr \vdots \cr y_n^{}(t) }, \quad
% {\bf f}(t,{\bf y}) =
% \pmatrix{f_1^{}(t,{\bf y})\cr \vdots \cr f_n^{}(t,{\bf y}) }. $$
% The dependent variables are thus $y_1^{},\dots,y_n^{}$, and
% we could write the system out explicitly like this:
% $$ y_1'(t) = f_1^{}(t, y_1^{}(t),\dots, y_n^{}(t)) , $$
% $$ \vdots $$
% $$ y_n'(t) = f_n^{}(t, y_1^{}(t),\dots, y_n^{}(t)) . $$
% For systems of small size we will usually find it
% convenient to use variable names without subscripts, such as $u,v,w$.
% The reason we can restrict the discussion to first order is that any
% higher-order ODE can be written as a system of first-order ODE\kern .5pt s by
% introducing additional variables, as outlined in
% Exercise 1.5.\ \ In {\sf FLASHI} nomenclature, any
% {\sf f\kern 1pt S$\,\cdots$} can be turned into
% an {\sf Fs$\,\cdots.$}\footnote{This doesn't
% mean that in practice, higher-order systems are necessarily reduced
% to first-order form.  For example, the Solar
% System Dynamics Group at the Jet Propulsion
% Laboratory in California solves systems of hundreds
% of ODE\kern .5pt s numerically to track the moon and
% planets and other bodies: when the media tell you there's going to
% be an eclipse,
% the computations probably came from JPL.\ \ The methods used
% by the SSD Group are based on a second-order formulation, without reduction
% to first order.}
% </latex>

%%
% <latex>
% To fully specify an ODE problem from (10.1), one normally
% needs $n$ additional conditions.  In an initial value problem,
% we prescribe the value of $\bf y$ at a fixed time such as $t=0$,
% $$ {\bf y}(0) = {\bf a} . \eqno (10.2) $$
% This amounts to specifying each of the $n$ components,
% $$ y_1^{}(0) = a_1^{}, ~\dots, \quad y_n^{}(0) = a_n^{} . $$
% In a boundary value problem, different components of $\bf y$ or
% combinations thereof may be
% specified at different points, for a total normally of $n$
% conditions all together.  Exactly what specifications may lead
% to a well-posed BVP for a system of equations
% is not always a straightforward matter.
% For IVPs, however, existence and uniqueness are guaranteed so
% long as $\bf f$ is continuous with respect to $t$ and $\bf y$ and Lipschitz
% continuous with respect to $\bf y$, as we shall prove in the next
% chapter (Theorem~11.2).
% </latex>

%%
% <latex>
% For example, consider the two-variable autonomous system
% $$ u' = -v - u(u^2+v^2), \quad v' = u - v(u^2+v^2). \eqno (10.3) $$
% Here are the trajectories of $u$ and $v$ for $t\in[\kern .3pt 0,30\kern .3pt]$
% with $u(0) = 1$, $v(0) = 0$.\footnote{For
% plots like this showing multiple components, we revert to
% Matlab's built-in color tables rather than the light/dark green and blue
% colors usually used in this book for linear/nonlinear IVPs and BVPs.}
% </latex>
ODEformats, N = chebop(0,30); N.lbc = [1; 0];
N.op = @(t,u,v) [diff(u)+v+u*(u^2+v^2); diff(v)-u+v*(u^2+v^2)];
[u,v] = N\0; plot([u v]), title('Fig.~10.1:~~Equation (10.3)',FS,11)
%%
% <latex> \vskip 1.02em </latex>
%%
% <latex>
% \noindent
% We see two interlocking oscillations, decaying rapidly in amplitude at
% first but then much more slowly.
% </latex>

%%
% <latex>
% Since (10.3) is autonomous, the whole future
% of the system is determined by the values of $u$ and $v$ at
% any given time.  Thus it makes sense
% to plot this solution in the {\bf phase plane},
% with $u$ on the horizontal axis and $v$ on the vertical axis.
% This is a straightforward generalization of
% the phase plane of the last chapter, where the two
% axes represented $y$ and $y'$.
% Here, however, trajectories need not point
% rightward in the upper half-plane.
% </latex>
arrowplot(u,v,CO,ivpnl,LW,.5,MS,4)
title('Fig.~10.2.~~Solution of (10.3) in the phase plane',FS,11)
axis equal, ylim([-.4 .6]), hold on, plot(0,0,'.k',MS,8)
xlabel('$u$',FS,10,IN,LT), ylabel('$v$',FS,10,IN,LT)
%%
% <latex> \vskip 1.02em </latex>
%%
% <latex>
% \noindent
% The dot at the center marks the fixed point of the
% system, corresponding to the constant solution
% $u(t) = v(t) = 0$.  In general, a {\bf fixed point} of
% an autonomous ODE ${\bf y}' = {\bf f}({\bf y})$
% is any point ${\bf y}_0^{}$ such that 
% ${\bf f}({\bf y}_0^{}) = {\bf 0}$, implying that
% the constant ${\bf y}(t) = {\bf y}_0^{}$ is a solution.
% The dynamics of this problem is angularly symmetric, as
% can be confirmed by converting $\bf f$ to polar coordinates.
% </latex>

%%
% As in the last chapter, we may think of the trajectory as
% following the arrows of a direction field:
quiver(N,[-.7 1.1 -.4 .6],'k','scale',2,'xpts',12,'ypts',8)
axis equal, ylim([-.4 .6]), hold off
title('Fig.~10.3.~~Same with quiver arrows added',FS,11)

%%
% <latex>
% \noindent
% One way to understand the motion is to 
% interpret the nonlinear terms of (10.3) as linear terms
% with varying coefficients.    For a positive constant
% $\kappa$, the system
% $$ u' = -v - \kappa u, \quad v' = u - \kappa v $$
% describes a simple oscillation with damping.  Locally, near
% any particular point $(u(t), v(t))$, (10.3) will behave
% in just this way with $\kappa = u(t)^2+v(t)^2$.  This
% explains the behavior of the two images above: we start
% with a damping coefficient $\kappa = 1$, but at the
% end of the trajectory, $\kappa$ is about 60 times smaller:
% </latex>
kappa = u(end)^2 + v(end)^2

%%
% <latex>
% We emphasized in the last chapter that phase plane or more
% generally phase space analysis (for $n\ge 3$)
% is appropriate only for a system that is autonomous,
% so that the state $\bf y$ at a fixed time
% determines fully the future solution ${\bf y}(t)$.\footnote{A
% nonautonomous system of size $n$ can, however, be written as an autonomous
% one of size $n+1$ by introducing a new dependent variable $\tau$
% equivalent to $t$, governed by the equation $\tau' = 1$ and
% the initial value $\tau(0) = 0$.  See Figure 9.9.}  A fundamental property 
% of phase space analysis is that 
% {\em distinct trajectories can never meet}
% (assuming $\bf f$ is continuous with respect to $t$ and $\bf y$ and
% Lipschitz continuous with respect to $\bf y$).
% This follows from the uniqueness of solutions, for
% if an autonomous ODE has two solutions ${\bf y}_1^{}(t)$ and ${\bf y}_2^{}(t)$
% such that ${\bf y}_1^{}(t_1^{}) = {\bf y}_2^{}(t_2^{})$
% for some $t_1^{}$ and $t_2^{}$, then if solutions are unique,
% we must have ${\bf y}(t_1^{}+t) = {\bf y}_2^{}(t_2^{}+t)$ for
% all $t$, positive and negative.  So the trajectories corresponding
% to $\bf y_1^{}$ and $\bf y_2^{}$ must be the same.
% </latex>

%%
% <latex>
% For another example, let us consider
% the two-variable system known as the {\bf predator--prey} or
% {\bf Lotka--Volterra} equations.
% We suppose an environment in which there is a prey species, say rabbits,
% whose population as a function of time is $u(t)$, and a predator
% species, say foxes, whose population as a function of time is
% $v(t)$.\footnote{Volterra's original application of 1926
% involved fish populations in the Adriatic
% Sea.  Danby writes in {\em Computer Modeling: From Sports
% to Spaceflight$\dots$ From Order to Chaos:} ``Presentations of
% this model include many pairs: birds and worms, wolves and deer,
% sharks and tuna, to name a few.  I have found that my students prefer
% foxes and rabbits, so I shall refer to these.  After all, who
% identifies with a worm?''}
% In the absence of foxes, the rabbit population grows
% exponentially thanks to abundant vegetation, but it is decreased by
% fatal encounters with foxes at a rate proportional to
% the product $uv$,
% $$ u'= u - uv. \eqno (10.4) $$
% The fox population, on the other hand, decays exponentially in the
% absence of rabbits to eat but expands when they are available,
% $$ v'= -\textstyle{1\over 5}v + uv. \eqno (10.5) $$
% (The constant $1/5$ is included to make the dynamics nontrivial.
% In a serious application
% one would have dimensional constants in front of each term, based
% on measurements of field data.)
% Combining the two equations with initial conditions, 
% we obtain the Lotka--Volterra IVP
% $$ u' = u-uv, ~ v'= -\textstyle{1\over 5}v + uv, \quad t\ge 0, ~ u(0) = u_0^{},
% ~v(0) = v_0^{}. \eqno (10.6) $$
% </latex>

%%
% <latex>
% \noindent
% A quiver plot looks like this:
% \noindent
% </latex>
N = chebop(0,60); N.lbc = [1;1];
N.op = @(t,u,v) [diff(u)-u+u*v; diff(v)+.2*v-u*v];
quiver(N,[0 1.5 0 3],'k','xpts',12,'ypts',8)
xlabel('rabbits, $u$',FS,10,IN,LT), ylabel('foxes, $v$',FS,10,IN,LT)
title('Fig.~10.4.~~Quiver plot for Lotka--Volterra equations (10.6)',FS,11)

%%
% <latex>
% This problem has periodic solutions,
% where the rabbit population oscillates between large and (very) small and
% the fox population oscillates between large and (not so) small,
% with a phase lag since a large rabbit population leads to a growing
% fox population.  Here is the periodic solution
% that evolves from the initial populations $u(0)=v(0)=1$.
% </latex>
[u,v] = N\0; plot([u v])
title(['Fig.~10.5.~~Lotka--Volterra ' ...
       '(= predator--prey) equations'],FS,11)
text(40,1.8,'foxes',IN,LT,FS,9)
text(36.5,0.7,'rabbits',IN,LT,FS,9)
hold on, plot([0 0;60 60],[.2 1; .2 1],'--k',LW,.6), hold off
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The period is about 17.5.
% </latex>
[~,m] = max(v,'local'); period = diff(m)

%%
% <latex>
% To understand the Lotka--Volterra dynamics, again it is helpful to
% interpret the nonlinear terms of (10.4)--(10.5) as linear terms
% with varying coefficients.   Equation (10.4) can be written
% $u' = (1-v) u$, showing that locally,
% $u$ will decrease exponentially for $v>1$
% (the fox curve lies above the upper dashed line)
% and increase exponentially for $v<1$ (it lies below
% the upper dashed line).
% Thus the rabbit population has an exponential
% growth tendency that shuts off as $v$ rises above~$1$.  
% Similarly, equation (10.5) can
% be written $v' = (u-1/5) v$, showing that the fox population has
% an exponential decay tendency but grows
% whenever there are enough rabbits, $u>1/5$ (the rabbit curve lies
% above the lower dashed line).
% </latex>

%%
% <latex>
% Here is the orbit just computed in the
% phase plane, winding around counterclockwise several times.
% </latex>
arrowplot(u,v,CO,ivpnl), hold on
quiver(N,[0 1.5 0 3],'k'), axis([0 1.5 0 3])
title('Fig.~10.6.~~Lotka--Volterra cycle in the phase plane',FS,11)
xlabel('rabbits, $u$',FS,10,IN,LT), ylabel('foxes, $v$',FS,10,IN,LT)
plot([.2 0;.2 1.5],[0 1; 3 1],'--k',LW,.6), hold off
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Let us superimpose several orbits in a single plot, corresponding to
% the solutions arising
% from various initial populations of rabbits.
% </latex>
for s = .4:.2:1.4
  N.lbc = [s; 1]; [u,v] = N\0;
  arrowplot(u,v,LW,.8,CO,ivpnl,MS,4,YS,.5), hold on
end
title('Fig.~10.7.~~Superposition of six orbits',FS,11)
axis([0 1.5 0 3])
xlabel('rabbits, $u$',FS,10,IN,LT), ylabel('foxes, $v$',FS,10,IN,LT)
plot([.2 0;.2 1.5],[0 1; 3 1],'--k',LW,.6), hold off
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The figures show that in the Lotka--Volterra model, the
% effects of the initial conditions last forever:
% some fox-rabbit oscillations
% have great amplitude while others move little.  At the
% center is obviously a fixed point, which we may confirm
% by setting $u'=v'=0$ in (10.6).  This gives the equations
% $$ u = uv, \quad \textstyle{1\over 5} v = uv, $$
% implying that these equations in fact
% have two fixed points, one at $(1/5,1)$ and the other at $(0,0)$.
% Let us mark them in the plot.
% </latex>
hold on, plot([0 .2],[0 1],'.k',MS,8), hold off
title('Fig.~10.8.~~Same with fixed points added',FS,11)
%%
% <latex> \vskip 1.02em </latex>


%%
% <latex>
% The Lotka--Volterra equations are perhaps the most famous
% system of two ODE\kern .5pt s, appearing in every ODE textbook. 
% Turning now to three dimensions, it is only
% fitting that we should begin with the most famous system of
% three ODE\kern .5pt s, the {\bf Lorenz equations.}
% In 3D, a new possibility opens up for reasons we shall
% explain in a moment: the dynamics of an ODE
% may be {\em chaotic}.
% The Lorenz equations are the archetypical chaotic system,
% consisting of three coupled first-order equations,
% $$ u' = 10(v-u), \quad v' = u(28-w)-v, \quad w' = uv -(8/3)w.
% \eqno (10.7) $$ 
% The coefficients $10$, $28$, and $8/3$ are arbitrary, and it is
% not necessary to use exactly these values.  However, so much
% study of the Lorenz equations has been based on
% these parameter choices that they have become standard.
% </latex>

%%
% <latex>
% Equations (10.7) may look complicated at first
% glance, but in fact, the nonlinearities here are rather
% simple: the equation for $u'$ is linear, and the
% equations for $v'$ and $w'$ are quadratic, their only
% nonlinearities consisting of the term $-uw$ in the equation for $v'$
% and the term $uv$ in the equation for $w'$.
% To make an IVP, we impose initial conditions such as
% $$ u(0)= v(0) = -15, ~~ w(0) = 20. \eqno (10.8) $$
% Here is the solution of (10.7)--(10.8) up to time
% $t=30$, showing the three components together.
% </latex>
N = chebop(0,30);
N.op = @(t,u,v,w) [diff(u)-10*(v-u); ...
    diff(v)-u*(28-w)+v; diff(w)-u*v+(8/3)*w];
N.lbc = [-15; -15; 20];
[u,v,w] = N\0; plot([u v w])
title('Fig.~10.9.~~Lorenz equations (10.7)--(10.8)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% To focus on just one component of the solution, here is a plot of $v(t)$.
% Note the characteristic property of chaos: apparent
% randomness in a deterministic system.  We shall discuss
% chaos more fully in Chapter~13.
% </latex>

plot(v,CO,ivpnl,LW,.7)
title('Fig.~10.10.~~Component $v(t)$',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% The Lorenz equations are autonomous, so they can be
% represented in a phase space, which will have three dimensions.
% Here is a phase space plot of the trajectory just computed.  
% The plot includes dots showing the fixed points, which are the
% solutions of the equations we obtain by setting
% $u'=v'=w'=0$ in (10.7),
% $$ v = u, \quad 28u - uw =v, \quad uv = 8w/3. $$
% There are three fixed points, one at
% $(u,v,w) = (0,0,0)$ and the others at
% $$ (u,v,w) = (\pm 6\sqrt 2, \kern 1pt \pm 6\sqrt2, \kern 1pt 27). $$
%%
% <latex> \vskip -2em </latex>
plot3(u,v,w,CO,ivpnl,LW,.4)
xlabel('$u$',FS,10,IN,LT), ylabel('$v$',FS,10,IN,LT)
zlabel('$w$',FS,10,IN,LT)
title('Fig.~10.11.~~3D phase space for the Lorenz eqs.',FS,11)
c = 6*sqrt(2);
hold on, plot3([0 c -c],[0 c -c],[0 27 27],'.k',MS,6)
view(-35,30), axis([-17 19 -25 25 0 45]), hold off
set(gcf,'position',[380 220 540 300])
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% Ideas related to phase spaces and their trajectories are
% at the heart of the rich subject of {\em dynamical
% systems,} and a great deal is known about the
% mathematics and application of these notions.
% The traditional method of this analysis is to begin
% by finding fixed points and then examine the behavior near
% these points by approximating the nonlinear problem locally 
% by a linear problem.  As we shall see in Chapter 14,
% linear analysis makes sense at any point of a phase
% space, not just a fixed point; locally, any smooth nonlinear
% ODE might as well be linear.   What makes the fixed points
% particularly interesting is that trajectories interact in
% nontrivial ways at such points, and to see the
% overall dynamics in a phase space, it is often fruitful to
% start from the fixed points and see how they may be connected
% by various trajectories.  See Chapter 15.
% </latex>

%%
% <latex>
% One of the fundamental things one can learn from this
% kind of geometrical reasoning is that trajectories can get tangled up in
% 3D, but not in 2D.\ \ In 2D, since trajectories cannot cross, tangling
% is not possible.  Therefore, chaos can appear in autonomous
% systems with $n\ge 3$ but never with $n=2$.\footnote{A
% precise formulation of this idea is provided by the
% {\em Poincar\'e--Bendixson theorem.}}  Another point
% is that since a linear system of ODE\kern .5pt s behaves in essentially
% the same way at every point of phase space, its trajectories look essentially
% the same everywhere.  So they can't get tangled up either, and
% chaos can never appear in a linear problem, no matter
% how many variables there are.  Since we have not
% given a mathematical definition of chaos, these observations
% must be regarded as heuristic, but they can be made precise.
% </latex>

%%
% <latex>
% Let us finish with an example of a three-variable problem
% that is not chaotic, analyzed in pp.~202--204 of
% Bender and Orszag, {\em Advanced Mathematical Methods for
% Scientists and Engineers,}
% $$ u' = vw, ~~ v' = -2uw, ~~ w' = uv. \eqno (10.9) $$
% It can be shown that along any trajectory, the quantity
% $u(t)^2 + v(t)^2 + w(t)^2$ is constant, and thus trajectories
% lie on spheres.  Here is a plot showing some of them.
% </latex>
close all, N = chebop(0,10);
N.op = @(t,u,v,w) [diff(u)-v*w; diff(v) + 2*u*w; diff(w)-u*v];
for theta = -1.5:.2:1.5
  N.lbc = [cos(theta); sin(theta); 0];
  [u,v,w] = N\0; plot3(u,v,w,CO,purple,LW,.5), hold on
  plot3(-u,v,w,LW,.5,CO,blue), plot3(w,v,u,CO,green,LW,.5)
  plot3(w,v,-u,'k',LW,.5)
end
plot3([0 0 0 0 1 -1],[0 0 1 -1 0 0],[1 -1 0 0 0 0],'.k',MS,8)
xlabel('$u$',FS,5,IN,LT), ylabel('$v$',FS,5,IN,LT)
zlabel('$w$',FS,5,IN,LT), set(gca,FS,6), axis equal, xlim([-2 2])
title('Fig.~10.12.~~Phase space for (10.9)',FS,11), hold off
set(gcf,'position',[380 220 540 300])
%%
% <latex> \vskip 1.02em </latex>
%%
% <latex>
% \noindent
% As explained by Bender and Orszag, this example is related to
% the mechanical phenomenon known as {\em Eulerian wobble} and can
% be used to explain why if you toss a book
% in the air, you can make it spin stably around
% the shortest axis or the longest, but not the middle one: it will
% somersault.  This is also called the {\em tennis racket theorem}
% or the {\em Dzhanibekov effect.}
% </latex>

%%
% <latex>
% \begin{center}
% \hrulefill\\[1pt]
% {\sc Application: SIR model for epidemics}\\[-3pt]
% \hrulefill
% \end{center}
% </latex>

%%
% <latex>
% A major advance in epidemiology dates to a 1927 study by W. O. Kermack and
% A. G. McKendrick
% of a 1906 outbreak of bubonic plague in Bombay, now Mumbai
% (``A contribution to the mathematical
% theory of epidemics,'' {\em Proceedings of the Royal Society of London A}.)
% Kermack and McKendrick proposed
% what is now called the {\em SIR model\/} for the spread of
% disease.  Although the SIR model
% is now considered too simplistic for quantitative
% tracking of real epidemics, it exhibits several important properties
% and remains the foundation for more comprehensive models.
% </latex>

%%
% <latex>
% This is a {\em compartment model,} in which one imagines
% members of a population making the transition continuously from a box
% labeled S for
% {\em susceptible} to another box~I for {\em infected} and eventually
% to a box R for {\em recovered} or {\em removed} (possibly by death).
% In nondimensional form the system is
% $$ S' = -\beta S I, \quad I' = (\kern .5pt
% \beta S-\gamma)I, \quad R' = \gamma I,
% \eqno (10.10) $$
% where $\beta,\gamma>0$ are parameters
% and $S(t)$, $I(t)$, and $R(t)$ denote the proportions of the population
% in the three compartments.
% These equations model the dynamics of the
% transitions ${\rm S}\to {\rm I}$ and ${\rm I}\to {\rm R}$.
% First, $S$ decreases and $I$ increases at the rate
% $\beta S I$ (so $S$ is monotonically nonincreasing).
% This product corresponds to infected people
% coming into contact with susceptible ones with
% infection rate~$\beta$ per contact (high for measles, low for AIDS).
% Second, $R$ increases and~$I$ decreases
% at the rate $\gamma I$ (so $R$ is monotonically 
% nondecreasing).  This represents the approximation that
% infected people come to the end of their infections, whether
% through death or recovery, in a random fashion at a steady rate (again high
% for measles, low for AIDS).
% </latex>

%%
% <latex>
% Note that (10.10) is nonlinear because of the product
% $SI$.  Note also how easily we can spot factors that
% are omitted from the model.  For example, (10.10) says nothing
% about the course of the infection; an infected patient is
% simply assumed to recover or die at random with a 
% probability that does not vary from day to day.  Also, the
% model assumes the populations are perfectly mixed,
% with no spatial or social aspect to the epidemic.  In actuality,
% spatial and socioeconomic factors can be very important,
% as was made famous by John Snow's study of the cholera outbreak
% linked to a certain public water pump in a crowded quarter
% of London in 1854.
% </latex>

%%
% The equations (10.10) have two obvious steady states.
% If $S=1$ and $I=R=0$, nobody is infected and the disease
% never gets going.  If $S=I=0$ and $R=1$, the whole
% population is recovered or removed and 
% the epidemic has passed into history.
% The nontrivial dynamics emerges when we begin with a nonzero
% infected fraction $I_0^{}$,
% $$ S(0)=1-I_0^{}, \quad I(0)=I_0^{}, \quad R(0)=0. $$ 
% Here we show the course of the epidemic with $\beta = 2$,
% $\gamma = 1$, and $I_0^{} = 0.0001$.
beta = 2; gamma = 1; I0 = 0.0001; N = chebop([0 30]); close all
N.op = @(t,S,I,R) [ diff(S) + beta*S*I
                    diff(I) - beta*S*I + gamma*I
                    diff(R) - gamma*I ];
N.lbc = [1-I0; I0; 0];
[S,I,R] = N\0;  plot([S,I,R]), hold on, plot(S+I+R,'--k'), hold off
xlabel('time $t$',IN,LT,FS,9), ylabel('proportion of population',FS,9)
title('Fig.~10.13.~~An SIR epidemic (10.10)',FS,11)
text(19.3,.07,'infected, $I(t)$',FS,9,IN,LT)
text(19.3,.27,'susceptible, $S(t)$',FS,9,IN,LT)
text(19.3,.86,'recovered or removed, $R(t)$',FS,9,IN,LT)
text(19.3,1.06,'total $S(t)+I(t)+R(t)$',FS,9,IN,LT)
%% 
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% The model may be simple, but there are
% some remarkable things to see in these curves.
% First of all, note how a disease
% may be active in a small fraction of a population for
% a long time, exponentially growing but at a rate so
% low that it may go unnoticed.\footnote{What do you
% suppose is happening in a carton of milk in your
% refrigerator when its ``best before'' date is still a week away?}
% </latex>

%%
% Why does it grow at all?  As mathematicians, we may say
% that the fixed point $(S,I,R) = (1,0,0)$ is
% unstable.  Looking at (10.10), we see that the equation
% giving exponential growth is the second one,
% $I' = (\beta S-\gamma) I$.  If $\beta S > \gamma$, the
% coefficient is positive and we can expect $I$ to grow exponentially
% at the rate $\beta S - \gamma$.  For our choice of parameters, this is indeed
% the case at $t=0$, with  $\beta S - \gamma = 1$.
% As epidemiologists, we can interpret this condition as follows:
% each newly infected person, on average, infects more than 1 additional
% person before ceasing to be infective.

%%
% Eventually, on the other hand, the epidemic shuts off.  As soon
% as $S$ falls to the point where $\beta S\le \gamma$, $I(t)$ begins
% to diminish and the epidemic starts fading out.  
% In this example this occurs at the point 
% $t_{\rm c}$ where $S(t_{\rm c}^{}) = 0.5$.  
tc = roots(S-0.5)

%%
% <latex>
% \noindent
% Note that although the infected population begins
% to diminish after this point, the epidemic has by
% no means yet run its course.  In the end, most of
% the population has had the disease:
% </latex>
Rinfty = R(end)

%%
% <latex>
% \noindent
% But not all.  This is a fascinating and important aspect
% of epidemiology.  For a population to be safe from a certain
% pathogen, it is not necessary that everybody be immune.  It suffices
% for a high enough fraction to be immune so that further infections
% will not spread.  This is the celebrated effect
% of {\em herd immunity.}  And here we realize that the two ``obvious''
% fixed points of this system
% noted earlier are not the only ones.  In fact,
% any values $S_\infty^{}$ and
% $R_\infty^{} = 1-S_\infty^{}$ in $[\kern .3pt 0,1]$ 
% may be combined
% with $I_\infty^{} = 0$ to make a fixed point.
% </latex>

%%
% As $\beta$ increases, the peak of $I(t)$ comes both sooner and
% higher:
for beta = [1.3 1.6 2 2.8]
    N.op = @(t,S,I,R) [ diff(S) + beta*S*I 
                        diff(I) - beta*S*I + gamma*I
                        diff(R) - gamma*I ];
    [S,I,R] = N\0; plot(I), hold on, [val,pos] = max(I);
    text(pos,val+.02,['$\beta =$ ' num2str(beta)],FS,9,IN,LT,HA,CT)
end
xlabel('time $t$',FS,9,IN,LT), ylabel('infected proportion',FS,9)
ylim([0 .35])
title('Fig.~10.14.~~SIR epidemic for different infection rates',FS,11)
%% 
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% Epidemiological equations are not just used to track biological diseases.
% Similar compartment models have been used to study the spread of
% computer malware as well as other items that travel across
% networks --- like a rumor, or a meme.  See our favorite reference
% at the end of Chapter~15.
% </latex>


%%
% <latex>
% \smallskip
% {\sc History.}  
% Systems of differential equations were studied by
% d'Alembert and Lagrange in the mid-18th century,
% but vector notation was not yet
% available in those days.  Vector formulations
% became customary with the work of Peano and others near
% the end of the 19th century.
% </latex>

%%
% <latex>
% \smallskip
% {\sc Our favorite reference.}  J. D. Murray's wide-ranging book
% {\em Mathematical Biology\/}
% did much to create the field of mathematical
% biology as we know it today.  In its later extension to two volumes
% (Springer, 2001), the book focuses in volume 1 on systems of
% ODE\kern .5pt s modeling biological systems, including insect
% populations, fishery management, tumor cell growth,
% divorce rates, biological oscillators,
% propagation of neural signals, and much more.  The biology is
% varied, and the dynamical effects in the ODE\kern .5pt s are equally varied, combining
% to give a deeply satisfying illustration of the power of
% mathematics in understanding our world.
% \smallskip
% </latex>

%%
% <latex>
% \begin{displaymath}
% \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl
% {\sc Summary of Chapter 10.}   Any system of ODE\kern .5pt s can be
% reduced to a first-order system of the form
% ${\bf y}'(t) = {\bf f}(t,{\bf y})$. 
% An $n$-variable autonomous system
% ${\bf y}'(t) = {\bf f}({\bf y})$ has an $n$-dimensional
% phase space.  If $n\ge 3$ and the equations are
% nonlinear, there is the possibility of chaos.
% \vspace{2pt}}}
% \end{displaymath}
% </latex>

%%
% <latex>\smallskip\small\parskip=1pt\parindent=0pt
% {\em \underline{Exercise $10.1$}.  Four bugs on a rectangle.}
% In Chapter 3, using the complex arithmetic trick, we
% considered a scalar ODE in which a lion chases an
% antelope, always moving directly towards it.
% Now consider a problem of four bugs starting at four
% corners of a rectangle, each chasing the one counterclockwise
% from it at speed $1$.
% {\em (a)} Suppose the bugs begin at the corners of
% the square $(\pm 1,\pm 1)$.
% Draw a figure showing the evolution to time $t=1.99995$, and
% report the locations at that time of the four bugs.
% (This part of the problem could be reduced
% to a scalar ODE --- indeed it can easily be
% solved analytically --- but the next part is necessarily a system.)
% {\em (\kern .7pt b)} Now suppose the bugs begin at the four corners
% of the rectangle $(\pm 1.01, \pm 0.99)$.  Answer the
% same questions as in part {\em (a)}.  To learn about
% the extraordinary paths followed by the bugs as $t\to 2$, see
% Chapman, Lottes, and Trefethen, ``Four bugs on a rectangle,''
% {\em Proceedings of the Royal Society,} 2011.
% \par
% {\em \underline{Exercise $10.2$}.  Oregonator.}
% The Oregonator, or Field--Noyes model, is an approximation to the
% dynamics of the Belousov--Zhabotinskii chemical reaction, which exhibits
% periodic oscillations of its reagents. (You can see a video of the reaction at
% {\tt https://youtu.be/wxSa9BMPwow?t=3m9s}.) In nondimensional
% form the equations are
% $$ \varepsilon u' = qv - uv + u(1-u),\quad
% \delta v' = -qv - uv + w, \quad w' = u-w, $$
% where $q$, $\varepsilon$, and $\delta$ are parameters.
% For this exercise use the values $q=8\times 10^{-4}$, $\varepsilon=0.05$,
% and $\delta=0.02$. (These are not physically realistic for the BZ
% reaction.  For realistic values, the system is considered ``stiff'' and
% requires a change in the numerical solver Chebfun uses, as illustrated
% in equation 72 of Appendix~B.)\ \ Solve
% the system with Chebfun for $u(0)=w(0)=0$, $v(0)=1$, and plot
% the three components as functions of time. (A particularly interesting way
% to see the dynamics is {\tt comet3(u,v,w)}.) 
% \par
% {\em Exercise $10.3$.  Nonautonomous.}
% Suppose that a first-order, two-variable ODE system has unique
% solutions for any initial data and one of its solutions is
% $u(t)= \sin(t),$ $v(t) = \sin(2t)$.  
% Sketch this orbit in the phase plane and explain
% why this ODE cannot be autonomous. 
% \par
% {\em \underline{Exercise $10.4$}.  Matrix eigenvalues via ODE isospectral flow.}
% Suppose ${\bf A} = [\, a~b\,;~b~c\, ]$ is a $2\times 2$ real symmetric matrix.
% The eigenvalues of ${\bf A}$ can be determined from
% $\lim_{t\to\infty}^{}a(t)$ and $\lim_{t\to\infty}^{}c(t)$ 
% in the ODE $a' = 2\kern .3pt b^2$, $b' = b(c-a)$, $c' = -2\kern .3pt b^2$.
% {\em (a)} Try this for the matrix $[\, 3~2\,;~2~1\, ]$, plotting
% deviations from their limiting values
% as a function of $t\in [\kern .3pt 0,2\kern .2pt ]$.
% How close are $a(t)$ and $c(t)$ at $t = 1$ and $2$ to their
% final values?
% {\em (\kern .7pt b)} Do the same for the 
% matrix $[1~2\,;~2~3]$.  What is qualitatively different in this case?
% {\em (c)} What are the fixed points of this ODE?\ \ (If
% you are curious to learn more about isospectral
% flows in linear algebra, including the case of $n\times n$ matrices
% and connections with standard algorithms for computing matrix
% eigenvalues,
% see Deift, Nanda, and Tomei, ``Ordinary differential equations and
% the symmetric eigenvalue problem,'' {\em SIAM Journal on Numerical
% Analysis,} 1983.)
% \par
% {\em Exercise $10.5$.  Energy conservation for Eulerian wobble.}
% Show that the quantity $u^2+v^2+w^2$ is independent of time for
% solutions of the system (10.9).
% \par
% {\em \underline{Exercise $10.6$}.  Four kinds of fixed points}
% (from Hairer, N\o rsett, and Wanner).
% {\em (a)} Draw a quiver plot of the system $u' = (u-v)(1-u-v)/3$,
% $v' = u(2-v)$.  Determine the four fixed points analytically and
% include them in the plot.
% {\em (\kern .7pt b)} Calculate the eigenvalues of the Jacobian at these points
% and show that one is a sink, one is a source, one is a saddle,
% and one is a spiral.
% \par
% {\em \underline{Exercise $10.7$}.  Phase plane for SIR model.}
% {\em (a)} Suppose a solution of (10.10) satisfies 
% $S(0)+I(0)+R(0)=1$.  Show that 
% $S(t)+I(t)+R(t) = 1$ for all $t$.
% Using this conservation property, we can eliminate $S$ or $I$ or
% $R$ from the system and thus obtain a two-variable system that can
% be analyzed in the phase plane.
% {\em (\kern .7pt b)} Elimination of $R$ is particularly trivial.
% Write down the resulting two-variable system.  Solve it
% for the same parameters
% as in Figure~10.13 and make an arrowplot of the trajectory in
% the $S$-$I$ plane.  Where are the fixed points in this plane?
% </latex>
