%% 3. First-order scalar nonlinear ODEs

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

%%
% The first thing to do with nonlinear problems is enjoy them.
% The chances are you can't solve a particular nonlinear
% ODE analytically, but there
% is a wonderful variety of effects to explore.
% Incidentally, a change of notation applies
% here: mathematicians usually write a linear differential operator
% without parentheses, $y \mapsto Ly$,
% but a nonlinear one with parentheses, $y \mapsto N(y)$.  The
% change of letters from $L$ to $N$ is also a reminder of nonlinearity.

%%
% For example, here is a basic linear IVP, really just an
% integral since $y$ appears only in the term $y'$:
% $$ y' = 3\kern -.9pt\cos(t), \quad y(0) = 0. \eqno (3.1) $$
% Its solution is the sine wave $y(t) = 3\kern -.9pt\sin (t)$.
ODEformats, N = chebop(0,20); N.lbc = 0; N.op = @(t,y) diff(y);
rhs = chebfun('3*cos(t)',[0,20]);
y = N\rhs; plot(y,CO,ivp), ylim([-4 4])
title('Fig.~3.1.~~Linear oscillation (3.1)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% What happens if we add a nonlinear function of $y$
% to the operator?  For example, let us construct an ODE
% that behaves like (3.1) when $y$ has a small amplitude
% but ``shuts off'' as $|y|$ increases.  An example
% of such an equation is
% $$ y' + |y|^2 y = 3\cos(t), \quad y(0) = 0.
% \eqno (3.2) $$
% (We write $|y|^2 y$ instead of $y^3$
% to suggest the idea of an amplitude-dependent coefficient $|y|^2$
% multiplying the usual $y$ term.)
% The next plot shows the new solution superimposed on the previous one.  Note
% that the curves are about the same at first, but diverge as the
% amplitudes grow larger.
% </latex>
N.op = @(t,y) diff(y) + y^3;
y2 = N\rhs; plot(y,CO,ivp), ylim([-4 4])
hold on, plot(y2,CO,ivpnl), hold off
title('Fig.~3.2.~~Oscillation with nonlinear damping (3.2)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% Here is a variant of the same nonlinear idea.  Suppose that instead
% of a ``penalty'' for large values of $|y|$ we impose a
% ``barrier,'' preventing $|y|$ from reaching the value $1$.
% Here is an equation with a logarithmic barrier:
% $$ y' - {1\over 2} \log(1-|y|) y = 3\cos(t).
% \eqno (3.3) $$
% The solution is smooth, though it doesn't look it (see Exercise 3.1),
% and confined to values $-1 < y < 1$.
N.op = @(t,y) diff(y) - 0.5*log(1-abs(y))*y;
y3 = N\rhs; plot(y,CO,ivp), ylim([-4 4])
hold on, plot(y3,CO,ivpnl), hold off
title(['Fig.~3.3.~~Nonlinear damping by a ' ...
       'logarithmic barrier (3.3)'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% Examples like these, scalar problems of first order, are
% rather limited.  The variety will grow in later chapters
% when we come to problems of higher order or with
% multiple variables.

%%
% <latex>
% Both (3.2) and (3.3) are written with nonlinearities
% of the form $g(y) y$ for some function $g$ (equal to $|y|^2$ and
% $-0.5 \log(1-|y|)$, respectively).  This form suggests
% the idea that locally, a nonlinear ODE should behave
% approximately linearly.  Near any time $t\approx t_0^{}$, it should
% be possible to approximate the solution by the solution of a linear equation
% derived, say, by series expansion.  This is indeed the case,
% at least if the coefficients are continuous,
% and the idea will be made precise in Chapter~14.
% </latex>
 
%%
% <latex>
% For a few special classes of first-order nonlinear ODE\kern .5pt s, 
% analytical solutions are available.
% The most important category of such problems are 
% the {\bf separable equations.}
% In the last chapter we solved $y' = a(t)y$ by writing it as
% $$ {dy\over y } = a(t) dt . $$
% If we generalize $y^{-1}$ to a continuous
% function $b(y)$ that does not change sign, we get the equation
% $$ b(y) dy = a(t) dt , $$
% and the two sides can be integrated as before to get a solution
% at least locally.  We record the result as a theorem.
% </latex>

%%
% <latex> 
% \smallskip
% {\em
% {\bf Theorem 3.1.  Solution of first-order
% separable scalar homogeneous IVP
% (\textsf{\textbf{FlaSHI}}).}
% Let $a(t)$ be a continuous function of\/ $t$ and
% let $b(y)$ be a continuous nonzero function of\/ $y$.
% A solution $y(t)$ of the problem
% $$ b(y) dy = a(t) dt , \quad y(0) = y_0^{} \eqno (3.4) $$
% satisfies the equation
% $$ \int_{y_0^{}}^y b(x) dx = \int_0^t a(s) ds. \eqno (3.5) $$
% }
% \smallskip   
% </latex>

%%
% <latex>
% Equation (3.4) has the form $y' = f(t,y)$ (of type {\sf FlaSH}),
% but not every equation $y' = f(t,y)$ can be separated in this way.
% However, another case where separation is possible is if
% $f$ is independent of $t$.
% </latex>

%%
% <latex> 
% \smallskip   
% {\em
% {\bf Theorem 3.2.  Solution of first-order
% autonomous separable scalar homogeneous IVP
% (\textsf{\textbf{FlASHI}}).}
% Let $b(y)$ be a continuous nonzero function of~$y$.
% A solution $y(t)$ of the problem
% $$ b(y) dy = dt , \quad y(0) = y_0^{} \eqno (3.6) $$
% satisfies the equation
% $$ \int_{y_0}^y b(x) dx = t. $$
% }
% \smallskip   
% </latex>

%%
% <latex>
% \noindent
% Similarly $y' = f(t,y)$ is separable if $f$ is independent
% of $y$ --- it is just an integral.
% </latex>

%%
% <latex>
% A prototypical example of a problem of the form (3.6) is the ODE
% $$ y' = y^\alpha, \eqno (3.7) $$
% where $\alpha$ is a constant.  Dividing by $y^\alpha$ gives
% $dy/y^\alpha = dt$ and hence
% $$ {y^{1-\alpha}\over 1-\alpha}  = t - t_{\rm b}^{} $$
% for some constant $t_{\rm b}^{}$ (the letter b stands for
% ``blowup'').  This simplifies to
% $$ y = \big[(1-\alpha)(t - t_{\rm b}^{})\big]^{1/(1-\alpha)}.  \eqno (3.8) $$
% For $\alpha>1$, this solution exhibits a phenomenon
% of mathematical and physical interest,
% {\em blowup in finite time}.
% For a specific illustration, consider the IVP
% $$ y' = y^2, \quad y(0) = 1. \eqno (3.9) $$
% From (3.8) with $\alpha=2$, we see that the solution is
% $$ y = {1\over 1-t},  \eqno (3.10) $$
% which diverges to $\infty$ at $t_{\rm b}^{} = 1$.
% Here is a plot of the solution up to the point where
% it reaches the value $100$.
% </latex>
N = chebop(0,1); N.op = @(t,y) diff(y) - y^2;
N.lbc = 1; N.maxnorm = 100;
y = N\0; plot(y,CO,ivpnl), axis([0 1.12 0 120])
hold on, plot([1 1],[0 120],'--k'), hold off
text(1.01,30,'$t_b = 1$','color','k',IN,LT)
title('Fig.~3.4.~~Blowup in finite time (3.9)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% Behaviorally, (3.9) and (3.10) are rich in potential applications
% and interpretations.
% The ODE (3.9) describes a process where $y$ increases not just
% in proportion to its current amplitude, but faster.  For example,
% one can imagine that $y(t)$ represents the temperature at time
% $t$ of a smoldering haystack that smolders faster as it gets
% hotter.  The singularity at $t=1$ corresponds to the hay 
% catching fire --- spontaneous combustion or ``thermal runaway''. 
% (This simple idea is refined in the Application of Chapter~18.)
% </latex>

%%
% Now let us vary the problem slightly and replace (3.9) by
% $$ y' = y + y^2, \quad y(0) = 1. \eqno (3.11) $$
% Since the values of $y'$ are now bigger,
% we expect that blowup will occur
% again and it will happen sooner.
% Again a solution can be obtained by separation of variables.
% We write
% $$ {dy\over y + y^2} = dt, $$
% and integration gives
% $$ \log\left( {y\over 1+y}\right) = t-t_{\rm b}^{}, $$
% that is,
% $$ {y\over 1+y} = e^{t-t_{\rm b}}, $$
% or equivalently
% $$ y = {1\over e^{t_{\rm b}-t} - 1}. \eqno (3.12) $$
% With the initial condition $y(0)=1$, we have $t_{\rm b}^{} = \log 2$,
% confirming that as predicted, the blowup is earlier than in (3.10).
N.op = @(t,y) diff(y) - y - y^2;
y = N\0; plot(y,CO,ivpnl), axis([0 1.12 0 120])
hold on, plot(log(2)*[1 1],[0 120],'--k'), hold off
text(.71,30,'$t_b = \log(2)$','color','k',IN,LT)
title('Fig.~3.5.~~Exponential blowup in finite time (3.11)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% We solved (3.11) by separation of variables.  
% As it happens, this is a special case of
% a more general class of ODE\kern .5pt s that can also be solved analytically.
% A {\bf Bernoulli equation}\footnote{Named after Jacob Bernoulli, the
% brother of Johann and uncle of Daniel.}
% is an ODE of the form
% $$ y' = a(t) y + b(t) y^p,  \eqno (3.13) $$
% where $p$ is a constant and $a(t)$ and $b(t)$ are
% given functions of $t$.  
% (If $p=0$, the problem is linear, and if $p=1$ it is
% linear and homogeneous.  We assume $p$ is not $0$ or $1$.)
% If $a$ or $b$ is nonconstant, then separation of variables
% will not work for (3.13), but there is a different approach that
% still succeeds.
% Let us multiply by $y^{-p}$ to get
% $$ y' y^{-p}  = a(t) y^{1-p} + b(t). $$
% If we now make the substitution
% $$ u = y^{1-p}, \quad u' = (1-p) y^{-p} y',   $$
% the equation becomes
% $(1-p)^{-1} u' = a(t) u + b(t)$, that is,
% $$ u' = (1-p) a(t) u + (1-p) b(t).  $$
% This is a linear ODE, which can accordingly be
% solved by an integrating factor as described in
% Theorem 2.3, or by the method of undetermined coefficients.
% </latex>

%%
% For example, let us generalize the blowup problems (3.9) and (3.11)
% a little further.  Consider the ODE
% $$ y' = y + t\kern .5pt y^2, \quad y(0) = y_0^{}. \eqno (3.14) $$
% Comparing with (3.11), we see that for $y_0^{}=1$, the amplification
% will be weaker here for $t\in [\kern .3pt 0,1)$, so we can expect blowup
% at a time $t_{\rm b}^{} > \log 2$.   
% Dividing by $y^2$ converts the equation to 
% $$ {y'\over y^2} = {1\over y} + t, \quad y(0) = y_0^{}, $$
% and the change of variables $u = y^{-1}$ converts this to
% $$ u' + u = - t, \quad u(0) = u_0^{} = {1\over y_0^{}}. $$
% Applying Theorem 2.3 or the method
% of undetermined coefficients, we find that the
% solution of this IVP is
% $$ u(t) = 1-t+ e^{-t} (u_0^{} - 1), $$
% that is,
% $$ y(t) = {1\over 1-t+ e^{-t} (y_0^{-1}  - 1)}. \eqno (3.15) $$

%%
% Examining (3.15), we see that $y_0^{}=1$, as in (3.9) and (3.11), is
% a special case in which the exponential term does not appear.
% For $y_0^{}=1$, therefore, the blowup occurs at exactly the same
% time $t_{\rm b}^{} = 1$ as in (3.10).  Larger $y_0^{}$ brings blowup sooner,
% and smaller $y_0^{}$ defers it to later.  Here are
% solutions up to $t=1$ or $y=100$, whichever comes first,
% for initial values $y_0^{} = 0.90,0.91, \dots, 1.00$.
N.op = @(t,y) diff(y) - y - t*y^2;
for y0 = 0.90:.01:1
  N.lbc = y0; y = N\0;
  plot(y,LW,.7,CO,ivpnl), axis([0.95 1 0 120]), hold on
  title(['Fig.~3.6.~~Bernoulli eq.\ (3.14), ' ...
         '$y(0) = 0.90,0.91,\dots,1.00$'],FS,11)
end
hold off
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% We have explored at some length the phenomenon of
% blowup of solutions to an ODE.\ \  Physically, blowup
% in finite time corresponds to an explosion or another
% feedback process running out
% of control.  Mathematically, it illustrates
% the phenomenon of {\bf nonexistence} of
% solutions to certain nonlinear ODE problems.
% Because of the blowup at $t_{\rm b}^{}=1$, for example, there
% exists no solution to the problem (3.9) on the
% interval $[0\kern .3pt, 2\kern .2pt ]$.\footnote{In
% Chebfun, if the computed solution
% hits the limit {\tt N.maxnorm}, then all further values are set to
% {\tt NaN} --- not-a-number.}
% There is a well-established general theory of existence and uniqueness
% of solutions to ODE IVPs, presented in Chapter 11.
% The fundamental result of this theory asserts that the IVP
% $$ y' = f(t,y), \quad t\in [\kern .3pt 0,\infty), \quad y(0) = y_0^{} $$
% is guaranteed to have a unique solution if $f$ is
% continuous with respect to $t$ and Lipschitz continuous with
% respect to $y$.\footnote{This means that there exists a constant
% $C$ such that for all $t$ and $y_1^{}, y_2^{}$ in the range
% of interest, $|f(t,y_2^{})-f(t,y_1^{})| \le C |y_2^{}-y_1{}|$.
% Note that (3.9) loses Lipschitz continuity
% as $|y|\to \infty,$ implying that existence of solutions
% to this problem can fail only if $|y|$ diverges.}
% </latex>

%%
% <latex>
% From the ODE (3.7) with $\alpha<1$, we can develop
% an example of {\bf nonuniqueness} that has a beautiful
% physical interpretation.
% Consider (3.7) with $\alpha = 1/2$,
% $$ y' = y^{1/2}, \quad y(0) = 0. \eqno (3.16) $$
% From (3.8) we obtain the solution
% $$ y(t) = {1\over 4} t^2. $$
% An equally valid solution, however, is $y(t) = 0$, and this is
% the one that Chebfun will compute.  Alternatively, a solution to
% (3.16) might ``get going'' at any time $t_0^{}\ge 0$:
% $$ y(t) = \cases{ 0 & $t\le t_0^{}$ 
% \cr\noalign{\vskip 4pt} {1\over 4}(t-t_0^{})^2 & $t\ge t_0^{}$} . $$
% Thus there are not just two possible solutions
% but an infinite family.  Here is a plot of four of them.
% </latex>
t = chebfun('t',[0 6]);
for t0 = 0:3
  t = chebfun('t',[0, t0+2.6]);
  y = 0.25*(t-t0)^2*(t>t0); plot(y,CO,ivpnl), hold on
end
axis([0 6 -1 3]), hold off
title('Fig.~3.7.~~Four solutions to (3.16)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% If $y\ge y_0^{}$ for some constant $y_0^{} > 0$,
% then the right-hand side of (3.16) is
% Lipschitz continuous with respect to $y$.  This implies that uniqueness
% can fail for this problem only at points with $y=0$, though for
% any value of $t$.
% In Chapter~11 we shall see examples where uniqueness fails
% at isolated points such as $y=t=0$.
% </latex>

%%
% <latex>
% We mentioned a physical interpretation, and this involves
% the {\em leaky bucket problem.}  Suppose a bucket of water has a hole in
% the bottom, so the water flows out.  After a certain time, all the water
% will be gone, and then the bucket remains empty for all time.
% Here is the nonuniqueness effect
% in the words of Corless and Jankowski:\footnote{See ``Variations on
% a theme of Euler,'' {\em SIAM Review,} 2016.  Nonuniqueness for the leaky
% bucket is also discussed, among other places, in volume 1 of Hubbard and
% West, {\em Differential Equations:\ A Dynamical Systems Approach.}}
% \begin{quotation}
% \noindent Given an empty bucket, there's no way to tell when it was full --- if
% it ever was.
% \end{quotation}
% The connection with (3.16) is provided by {\em Torricelli's Law} of 1643 (Exercise
% 3.12).
% If $y>0 $ is the height of water in a leaky bucket, then
% $y$ decreases at a rate governed by the equation
% $$ y' = - C y^{1/2} $$
% for an appropriate constant $C$.  This means that if we take $t$ to be
% time measured {\em backward} from the present, the equation is
% $$ y' = C y^{1/2},$$
% which is essentially (3.16).
% So Figure 3.7, with time reversed, can be interpreted
% as a picture of the leaky bucket.
% </latex>

%%
% <latex>
% The last chapter was linear and this one is nonlinear, but both
% have been restricted to scalar problems of first order.
% Before turning to higher order problems and systems of equations,
% we want to illustrate an invaluable trick for simplifying
% computations involving particles in a plane: the use of
% complex arithmetic.  A particle moving in the $x$-$y$ plane
% has two coordinates, $x(t)$ and $y(t)$.  This suggests an
% ODE with two dependent variables, but if we define
% $z(t) = x(t) + iy(t)$, we have just a single, scalar variable
% $z(t)$.  This can be remarkably convenient.   One reason is that
% many particle interactions depend on distances
% between points in the plane, and a distance in the $x$-$y$ plane
% can be regarded as the absolute value of a complex vector.
% </latex>

%%
% Our Application illustrates this use of a complex variable.

%%
% <latex>
% \begin{center}
% \hrulefill\\[1pt]
% {\sc Application: classic pursuit problems}\\[-3pt]
% \hrulefill
% \end{center}
% </latex>

%%
% <latex>
% Suppose an antelope runs with speed 1 along the vertical line $x=1$,
% starting at $(1,0)$ and going in the positive $y$-direction. A lioness
% starts at $(0,0)$ and pursues the antelope, always moving directly toward
% the antelope at a fixed speed $C$.\footnote{The cast
% of characters is variable, as described in
% P. J. Nahin, {\em Chases and Escapes: The Mathematics
% of Pursuit and Evasion,} Princeton, 2006.
% One may have a dog chasing a rabbit, a hawk chasing a sparrow, or, in
% the original treatment by Pierre Bouguer in 1732, a pirate ship
% chasing a merchant ship.}\ \ What path does the lioness follow,
% and when and where does she catch the antelope?
% </latex>

%%
% <latex>
% The antelope's path is given as the function
% $a(t) = 1 + i\kern .3pt t\kern .5pt$; the 
% unknown in this problem is the path of the lioness. 
% In complex arithmetic, we can regard this as a
% function $z(t)$ given by the IVP
% $$ z' = C {a(t)-z(t)\over |a(t) - z(t)| }, \quad
% z(0) = 0. \eqno (3.17) $$
% Here is a plot showing the lioness's track up to time
% $t=4$ for $C = 0.5$, $1$, and~$1.1$.
% </latex>
tmax = 4; a = chebfun('1+1i*t',[0 tmax]);
N = chebop(0,tmax); N.lbc = 0; CC = [.5 1 1.1];
for j = 1:3
  C = CC(j); N.op = @(t,z) diff(z) - C*(a(t)-z)/abs(a(t)-z);
  subplot(1,3,j), plot(a,':r'), hold on, plot(a(end),'.r',MS,10)
  z = N\0; arrowplot(z,CO,ivpnl,YS,0.3), hold off
  text(0.75,4.6,['$C = ' num2str(C) '$'],FS,9,HA,CT)
  axis([0 1.5 0 5]), grid on
end
title('Fig.~3.8.~~The lion chases the antelope\kern .25in',FS,11,HA,RT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% With $C=0.5$, the lioness will obviously never catch
% the antelope, which soon sprints out of reach. 
% With $C=1$, she still never quite makes the catch.
% The separation distance is
% $0.50000015$ at $t=4$, converging exponentially
% to $1/2$ as $t\to\infty$.  For any $C>1$, on the other hand,
% the catch will take place.  We have stopped this experiment 
% before that point since there is a singularity involved.
% </latex>

%%
% <latex>
% Of course, the antelope may zig-zag.  Here is another
% run with $C=0.5$ and~$1$
% in which she makes a $90^\circ$ right turn at $t=2$.
% As $t\to\infty$ in the latter case, the separation distance 
% approaches $0.2653\dots .$
% </latex>
a = chebfun({'1+1i*t','-1+2i+t'},[0 2 4]);
N = chebop(0,tmax); N.lbc = 0; CC = [.5 1]; clf
for j = 1:2
  C = CC(j); N.op = @(t,z) diff(z) - C*(a(t)-z)/abs(a(t)-z);
  subplot(1,2,j), plot(a,':r'), hold on, plot(a(end),'.r',MS,10)
  z = N\0; arrowplot(z,CO,ivpnl,YS,1.4), hold off
  text(1.85,2.7,['$C = ' num2str(C) '$'],FS,9,HA,CT)
  axis([0 3.5 0 3]), set(gca,'ytick',0:3), grid on
end
title('Fig.~3.9.~~The antelope takes evasive action\kern -.1in',FS,11,HA,RT)
%%
% <latex> \vskip 1.02em </latex>

%%
% This pursuit problem has only a single unknown trajectory,
% making it a scalar problem in complex arithmetic and
% thus fitting the theme of this chapter.
% The main use of complex arithmetic, however, is for tracking systems
% of multiple particles in the plane, and we shall use this method
% to compute planar orbits of planets in Chapter 13 and electrons
% in Exercises 13.5 and 19.4,
% as well as looking at a multi-particle pursuit problem in
% Exercise 10.1.


%%
% <latex>
% \smallskip
% {\sc History.}  Nonlinear equations have been
% part of the study of ODE\kern .5pt s from the beginning.
% For the first thirty years or so, until about
% 1700, the equations were of first order.  Then
% higher order equations joined the discussion.
% </latex>

%%
% <latex>
% \smallskip
% {\sc Our favorite reference.}  Even if you don't read German,
% the book {\em Differentialgleichungen:\ L\"osungsmethoden
% und L\"osungen} by Erich Kamke
% (Springer Fachmedien Wiesbaden, 1977) is extraordinary.
% This book, which appeared in many editions starting in 1942,
% features a collection of 1600 numbered
% examples of ODE\kern .5pt s with their solutions --- the first 576 of
% them, filling 103 pages, corresponding to first-order equations.
% It is a monument to the
% knowledge of ODE\kern .5pt s accumulated during their first 300 years
% and a poignant indication of how different the world was
% before computers.
% </latex>

%%
% <latex>
% \def\equals{\kern 1pt =\kern 1.5pt}
% \begin{displaymath}
% \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl
% {\sc Summary of Chapter 3.}  
% Some first-order nonlinear ODE\kern .5pt s can be solved by separation of
% variables, and Bernoulli equations, of the form
% $y' = a(t)y + b(t)y^p$, can be solved by the change of variables
% $u = y^{1-p}$.  Most other first-order
% nonlinear problems cannot be solved
% analytically.  The ODE\/ $y' = f(t,y)$ has a unique solution 
% for all\/ $t$ if\/ $f$ is continuous with respect to\/ $t$ and Lipschitz
% continuous with respect to\/ $y$.  If these conditions do not hold,
% existence and/or uniqueness may fail. \vspace{2pt}}}
% \end{displaymath}
% </latex>

%%
% <latex>
% \small\parindent=0pt\parskip=1pt
% {\em \underline{Exercise $3.1$}.  Smoothness of the solution to the barrier problem.} 
% Despite appearances in Fig.~3.3, the solution $y(t)$ of (3.3) is smooth.
% Confirm this by plotting $y,$ $y',$ and $y''$. 
% What are the minimum and maximum values of these three
% functions over the interval
% $[1,20\kern.3pt]$ (to exclude initial transients)? 
% (You can use a command like \verb|y{1,20}| to restrict a chebfun
% to a subinterval.)
% \par
% {\em Exercise $3.2$.  Analytical solutions via clever substitutions.} 
% One method used by experts in analytical solution of ODE\kern .5pt s is
% to change variables.  Find the general solutions of
% the following problems
% analytically using the substitutions indicated:
% {\em (a)} $y' = e^{t-y}-e^{\kern .3pt t}$ ($u = e^y$),
% {\em (\kern .7pt b)} $t\kern .3pt y' = y(\log(t\kern .3pt y)-1)$ ($u = t\kern .3pt y$),
% {\em (c)} $2t\kern .3pt yy' = y^2-t$ ($u = y^2$).
% \par
% {\em Exercise $3.3$.  A solution with compact support.} 
% Consider the IVP $y' = -t y/|y|^{1/2},$ $y(0) = 1$, where
% $|y|^{1/2}$ represents the positive branch of the square root.
% {\em (a)} Find analytically the (unique) solution for $t\in [\kern .3pt 0,1]$.
% {\em (\kern .7pt b)} Find analytically the (unique) solution for $t\in [\kern .3pt 0,\infty)$.
% \par
% {\em \underline{Exercise $3.4$}.  Spherical flame in microgravity.}
% (Adapted from section 7.9 of C. B. Moler, {\em Numerical Computing with
% MATLAB,} SIAM, 2008.
% A video of a growing flame sphere in microgravity can be found at
% {\tt https://goo.gl/nQ5Vxd}.)
% In the absence of gravity, a flame takes a nearly spherical shape.
% Oxygen, which fuels the flame, enters the sphere at a rate proportional
% to its surface area. Combustion consumes the oxygen at a rate
% proportional to the volume. In appropriate units the radius
% $r(t)$ is approximately governed by the ODE
% $dr/dt = r^2-r^3$.  We assume that $r(0)=r_0^{}$.
% {\em (a)} Show that the solution satisfies the implicit equation
% $\log(r/(1-r)) - 1/r = t + C$,
% and give a formula for $C$ in terms of the initial condition.
% {\em (\kern .7pt b)} Show that the time $t_h^{}$ at which
% $r$ takes the value $1/2$ is approximately $1/r_0^{}$ as
% $r_0^{}\to 0$. (Note: $1/x$ dominates $\log x$ as $x\to 0$.)
% {\em (c)} On one graph, plot numerically obtained solutions for $0\le t\le
% 1500$ and $r_0^{}=10^{-1},10^{-2},10^{-3}$.  As in Exercise 3.1, plot
% a zoom of the
% figure for $r_0^{} = 10^{-3}$ to confirm that despite appearances, this
% solution is smooth.
% \par
% {\em Exercise $3.5$.  Multiple routes to the same solution.} 
% Consider the ODE $y' + \sin(y) = 0$. {\em (a)} Find the general
% solution analytically using separation of variables.
% {\em (\kern .7pt b)} Find it again by interchanging the independent and
% dependent variables as in Exercise 2.3.
% \par
% {\em Exercise $3.6$.  Some nonlinear problems.} 
% Suppose $y(0)=1$.  Determine $y(1)$ analytically if
% {\em (a)} $y' = y^{3/2}e^{\kern .3pt t}$, {\em (\kern .7pt b)} $(t+1)y' +3y = 0$, {\em (c)} $yy' = t$.
% \par
% {\em Exercise $3.7$.  Fixed points and stability.} 
% A number $y_*^{}$ is a {\bf fixed point} of an autonomous ODE 
% $y' = f(y)$ if $f(y_*^{}) = 0\kern .7pt ;$
% it is {\bf stable} if $f'(y_*^{})<0$
% and {\bf unstable} if $f'(y_*^{})>0$.  Find all fixed points and determine
% their stability or instability for 
% {\em (a)} $y' = y + y^2$ (eq.~(3.11)),
% {\em (\kern .7pt b)} $y' = y^2 - 1$,
% {\em (c)} $y' = y - y^2$ (the logistic equation; see also Exercises
% 3.15 and 3.16), and {\em (d)} $y' = \sin(y)$.
% \par
% {\em \underline{Exercise $3.8$}.  Stable and unstable fixed points.} 
% {\em (a)} Now explore equation {\em (c)} of the last
% exercise on the computer by making a plot of the trajectories emanating
% from $y(0) = -1,-0.8, \dots, 1.8, 2$ on the interval
% $t\in [\kern .3pt0,4]$.  You will need to use the {\tt N.maxnorm}
% feature since otherwise some curves will blow up.
% {\em (\kern .7pt b)} Similarly, explore equation {\em (d)} of the last
% exercise with a plot of the trajectories emanating
% from $y(0) = -15,-14,\dots, 15$ on the interval
% $t\in [\kern .3pt0,4]$.
% \par
% {\em \underline{Exercise $3.9$}.  Ghost of a fixed point.} 
% Plot the solution of $y' = 1 - a \sin(y)$, $y(0)=0$ for
% $t\in [\kern .3pt 0,200\kern.3pt]$ with $a=0.9$, $0.99$, and $0.999$.
% In each case use {\tt roots(y-2*pi)} to determine the value
% of $t$ for which $y(t) = 2\pi$.
% \par
% {\em \underline{Exercise $3.10$}.  The Lambert W function.} 
% {\em (a)}
% The {\em Lambert W function} is a function $W(t)$ defined
% by the functional equation $W(t\kern .3pt e^{\kern .3pt t}\kern .2pt ) = t$.  It satisfies the
% differential equation $W'(t) = W(t)/(t+t\kern .3pt W(t))$ and takes the
% value $W(e) = 1$, where $e = 2.718\dots.$  Use this information
% to solve an ODE to compute the value $W(1)$.
% {\em (\kern .7pt b)} The number just computed probably matches the
% exact value of $W(1)$ (readily found on the Web) to about
% 10 digits of accuracy, because that is Chebfun's default accuracy
% for ODE solutions.  To get more digits, execute
% \verb|cheboppref.setDefaults('ivpAbsTol',1e-14)| before making
% the calculation (try {\tt chebfunpref} and {\tt cheboppref} and
% see Chapter 8 of the {\em Chebfun Guide} for more information).
% How accurate is the new value?
% Afterwards, return to the usual defaults by executing
% \verb|cheboppref.setDefaults('factory')|.
% \par
% {\em \underline{Exercise $3.11$}.  Antelope on a circle or square.} 
% As in Figs.~3.8 and 3.9, suppose a lion begins
% at $z=0$ at $t=0$ and chases the antelope with
% equal speeds, i.e., $C=1$ in (3.17).
% {\em (a)} Suppose the antelope runs around the unit circle with
% position $a(t) = e^{it}$.  Plot the lion's trajectory over the
% interval $t\in [\kern .3pt 0, 2\pi]$.
% Plot the distance between the two animals as a function of $t$.  At roughly
% what time $t$ will the distance fall below $0.01\kern .5pt?$
% {\em (\kern .7pt b)} Suppose the antelope runs around the unit square, which
% you can construct with
% \verb|a = 1 + cumsum(round(exp(pi*.25i*(t+2))/sqrt(2)))| if 
% {\tt t} is a chebfun for $t$ defined over the time interval of interest.
% Plot the lion's trajectory over the
% interval $t\in [\kern .3pt 0, 8]$.  
% Again plot the distance as a function of $t$ and estimate when this
% will fall below $0.01$.
% \par
% {\em Exercise $3.12$.  The leaky bucket problem ---
% Torricelli's law (1643).}
% A cylindrical tank has a hole at the bottom.  Water flows
% out, making the height $y(t)$ of the water decrease from its
% initial value $y_0^{}$ to $0$ after a certain time.  Derive the ODE for
% this process by considering energy, as follows.  The water in
% the tank has a certain {\em potential energy} determined by $y(t)$
% and hence decreasing at a rate determined by $y$ and $y'$.
% As water leaves the
% hole, this is converted to an equal amount of {\em kinetic energy} 
% determined by $y'$ and hence increasing at a rate determined by
% $y'$.  By balancing these two, explain why the ODE has the
% form $y' = -Cy^{1/2}$.  Solve the equation analytically
% and show that $y(t) = 0$ is reached at a finite time $t$.
% If $t_0^{}$ is doubled, what effect will this have on the drainage time?
% \par
% {\em \underline{Exercise $3.13$}.  Blowup equation with a
% complex coefficient.}
% Consider equation (3.9) except with a complex coefficient:
% $y' = Cy^2$, $y(0) = 1$, where $C$ is a constant with a nonzero
% imaginary part.
% {\em (a)} Write down the analytical solution, valid for all $t$.
% {\em (\kern .7pt b)} Plot the solutions corresponding to $C=1+1i$ and $C = 1+0.1i$.
% \par
% {\em \underline{Exercise $3.14$}.  Complex nonlinear oscillator.}
% The equation
% $y' = iy + 0.1(1-|y|^2) y$ might be regarded as a kind
% of complex, first-order analogue of the van der Pol equation.
% {\em (a)} Compute and plot the solution for $y(0)=0.1$ and
% $t\in [\kern .3pt 0,100\kern .5pt]$.
% {\em (\kern .7pt b)} Use $\tt roots(real(y))$ to determine the period
% of the oscillation, and explain why this is the value that appears.
% \par
% {\em \underline{Exercise $3.15$}.  Logistic equation.}
% Positive solutions of $y' = y$ grow exponentially, but 
% positive solutions of
% $y' = (1-y/Y)y$, where $Y>0$ is a constant known as
% the {\em carrying capacity}, asymptote to $y = Y$ as $t\to\infty$.  
% {\em (a)} Solve the equation analytically for $y(0) = y_0^{}$.
% {\em (\kern .7pt b)} Solve it numerically with $Y=5$ for $t\in [\kern .3pt 0,5]$
% and make a plot of trajectories with
% $y_0^{} = -5,-4,\dots, 10$.
% (This equation was derived by Pierre-Fran\c cois Verhulst in the
% 1830s as a model of population growth after he read Malthus's
% {\em Essay on the Principle of Population.})
% \par
% {\em \underline{Exercise $3.16$}.  Logistic equation with harvesting.}
% Suppose a population is naturally governed by logistic growth but is regularly
% harvested, obeying $y'=y(1-y/Y)-H$ for positive constants
% $Y$ and $H$. 
% {\em (a)} Show theoretically that the existence of a steady real solution
% satisfying $y'=0$ requires $H<Y/4$. 
% (b) Using Chebfun, plot
% solutions over $[\kern .3pt 0,11]$ for $y_0^{}=5$, $Y=5$, and
% $H=0.75,1,1.25,1.5$. How do you interpret the solution with $H=1.5\kern .5pt$?
% \par
% {\em \underline{Exercise $3.17$}.  Blowup equation with a $t$-dependent correction.}
% Positive solutions of $y' = y^2$ blow up to $\infty$ in finite
% time, but $y' = y^2 - t^2$ has different behavior.
% Plot solutions for $t\in[-2,2\kern .3pt]$ with $y(-2) = y_0^{} = 
% -2, -1.9, -1.8, \dots, 2$.
% (Use {\tt N.maxnorm = 3}.) 
% Describe in words what trajectories appear to look like as $t\to\infty$
% as a function of $y_0^{} \in [-2,2\kern .3pt]$.  Approximately what is
% the critical value of $y_0^{}$ in this interval?
% </latex>
