%% 16. Multiple solutions of nonlinear BVPs
% <latex>
% \setcounter{page}{196}
% </latex>

%%
% In the subject of linear algebra, existence and uniqueness
% are straightforward.  A scalar linear equation $ay = b$ has a unique
% solution for any $b$ if $a\ne 0$, whereas if $a=0$, it either
% has infinitely many solutions (if $b=0$) or 
% no solution at all (if $b\ne 0$).
% For a system of $n$ linear equations ${\bf A}{\bf y}={\bf b}$, where ${\bf A}$
% is an $n\times n$ matrix, the situation is a generalization
% of the same alternative.  If ${\bf A}$ is 
% nonsingular, there is a unique solution for any $\bf b$, whereas
% if ${\bf A}$ is singular, there are infinitely many solutions if
% ${\bf b}\in \hbox{range}({\bf A})$ and 
% no solutions at all if ${\bf b}\not\in \hbox{range}({\bf A})$.

%%
% Nonlinear algebraic problems, by contrast, can do almost anything.
% Consider first a nonlinear equation involving a scalar real variable $y$,
% which without loss of generality we can write with a
% zero right-hand side as
% $f(y) = 0$.  If $f(y) = \exp(y)$ there
% is no solution; if $f(y) = y+\exp(y)$ there is one solution; and
% if $f(y) = \hbox{round}(y)$ there
% are infinitely many solutions, namely all the numbers
% in the interval $(-0.5,0.5)$, together with perhaps $-0.5$ or $0.5$ depending
% on exactly how you define the ``round'' function.
% More important, it may happen that there are
% multiple solutions that are separated from one another.
% These may be finite in number, as with $f(y) = y^2-1$, whose solutions
% are $y=\pm 1$, or they may be
% infinite in number, as with $f(y) = \sin(\pi y),$ whose solutions
% are all the integers. The function $f(y) = \sin(\pi y) + 0.01y^2$
% illustrates another possibility, that a problem
% may have quite a few isolated solutions without having
% infinitely many.
ODEformats
f = chebfun(@(y) sin(pi*y)+0.01*y.^2,[-15 15]); plot(f,'k')
hold on, r = roots(f); plot(r,f(r),'.r',MS,12), hold off
grid on, ylim([-2 4])
title(['Fig.~16.1.~~A nonlinear function $f(y)$ ' ...
       ' with many zeros'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% All this is for a scalar equation $f(y)=0$.  For a nonlinear system
% of equations ${\bf f}({\bf y}) = {\bf 0}$, further possibilities arise.
% </latex>

%%
% Now let us look at differential as opposed to algebraic equations.
% What can we say about the possibility of nonuniqueness,
% that is, of problems with multiple solutions?

%%
% <latex>
% For linear ODE\kern .5pt s, we have already made the main observations.
% Theorems 2.1--2.4 asserted uniqueness of solutions
% to first-order scalar linear IVPs, and these results carry over to
% IVP linear systems and equations of
% higher order.  For linear BVPs, we saw in Chapter 6
% that the key matter is whether or not a problem has
% an {\em eigenfunction}.
% If not, there is a unique solution.  If so, there are either
% no solutions or a continuum of solutions.\footnote{For example,
% $y''+y = 0$ with boundary condition $y(0)=0$ has the general solution 
% $y(x) = A\sin(x)$.  On any interval $[\kern .3pt 0,L]$ where $L$
% is not an integer multiple of $\pi$, there is a unique solution for
% any $b$ when
% the second boundary condition $y(L) = b$ is specified.
% If $L$ is an integer multiple of $\pi$, on the other hand,
% there are infinitely many solutions if $b=0$ and no solutions
% if $b\ne 0$.}
% </latex>

%%
% <latex>
% For nonlinear ODE IVPs, uniqueness is again usually not
% a problem.  Theorem~11.1 asserted existence and
% uniqueness for any problem $y'(t) = f(t,y)$
% if~$f$ is continuous with respect to~$t$ and Lipschitz
% continuous with respect to~$y$, and Theorem 11.2 made the
% analogous statement for systems of IVPs.\ \ For uniqueness to fail
% for an IVP, $f$ must lack these continuity properties,
% as in the example $y'= y^{1/2}$ of equation (3.16) or the
% further examples discussed in Chapter~11.
% </latex>

%%
% <latex>
% When it comes to nonlinear ODE BVPs, however, anything
% is possible. We saw illustrations of nonlinear nonuniqueness
% with equation (5.10) and Exercise~9.2, and the BVPs associated
% with the nonlinear pendulum of Chapter 9 also have
% nonunique solutions, though we did not mention
% that there.  Phenomena of this kind
% are the subject of this and the next two chapters.
% </latex>

%%
% To make a start, consider the linear BVP
% $$ y'' = - y, \quad x\in [\kern.3pt 0,1], ~ y(0) = y(1) = 1. 
% \eqno (16.1) $$
% This problem has no eigenfunction, so there is a unique solution.
N = chebop(0,1); N.lbc = 1; N.op = @(x,y) diff(y,2) + y;
N.rbc = 1; y = N\0; plot(y,CO,bvp), ylim([0.9 1.3])
title(['Fig.~16.2.~~Linear BVP (16.1):  maximum $' ...
       num2str(max(y)) '$'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Now consider the nonlinear variant in which the right-hand
% side of (16.1) is replaced by its cube,
% $$ y'' = - y^3, \quad x\in [\kern.3pt 0,1], ~ y(0) = y(1) = 1. 
% \eqno (16.2) $$
% (A quintic as opposed to cubic ``nonlinear spring law'' appeared 
% in eq.~(4.8).)
% If we give this problem to Chebfun, a solution is produced that
% looks approximately like the last one, but about 10\% larger.
% </latex>
N.op = @(x,y) diff(y,2) + y^3;
y1 = N\0; plot(y1,CO,bvpnl), ylim([0.9 1.3])
title(['Fig.~16.3.~~Nonlinear BVP (16.2):  maximum $' ...
     num2str(max(y1)) '$'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% However, this is not the only solution to (16.2).  
% First let us present another solution; then we shall explain
% how we got it.
% </latex>
x = chebfun('x',[0 1]); N.init = 1-25*(x-x^2);
y2 = N\0; plot(y2,CO,bvpnl), ylim([-6 4])
title('Fig.~16.4.~~Another solution to (16.2)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% The last two figures show that there exist
% at least two distinct solutions
% to (16.2).  Now this is not a book of numerical analysis, but
% a word must be said at this point about algorithms.  In every area of
% computational mathematics, to solve nonlinear problems, it is usually
% necessary to use some kind of iteration, in which the
% solution is approached via a sequence of linear problems.
% The prototypical iteration
% is Newton's method, and Chebfun uses a version
% of Newton's method to solve nonlinear BVPs.\footnote{When
% applied like this to find functions as opposed to just numbers,
% Newton iteration is also called {\em Newton--Kantorovich iteration.}}
% </latex>

%%
% <latex>
% If a problem has more than one solution, which one will
% an iteration like Newton's method converge to?  The first
% thing to be said is that sometimes, it may not converge 
% at all, and there is a large subject of
% nonlinear numerical optimization that aims to improve matters
% in this respect.  When it does converge, however, the
% solution it converges to is often one that is close to the
% {\em initial guess} employed by the iteration.
% Every Newton iteration starts from some initial guess or
% other, even if it is the zero function.  In the case of
% Chebfun, if the user does not specify an initial guess
% explicitly, then the iteration starts from a simple polynomial
% in the variable $x$ constructed to match the boundary conditions,
% and this is what Chebfun did to obtain the solution
% plotted in Figure~16.3.
% </latex>

%%
% <latex>
% To get the second solution of (16.2), we overrode the 
% default by specifying a different initial
% guess in the field \verb|N.init|.
% This function, $1-25(x-x^2)$, was chosen because it has
% approximately the right shape.
% There are few guarantees in this business, but Chebfun
% duly converged to the solution plotted in Figure~16.4.
% Here we superimpose the initial guess on the plot.
% </latex>
hold on, plot(N.init,'--r'), hold off
title('Fig.~16.5.~~Solution shown with initial guess',FS,11)
%%
% <latex> \vskip 1.02em </latex>


%%
% <latex>
% As an autonomous scalar equation of second order, (16.2) can
% be examined in the $y$-$y'$ phase plane.  Here is an image showing the
% phase plane trajectories of the two solutions just computed.
% </latex>

y1p = diff(y1); y2p = diff(y2);
plot(y1([0 1]),y1p([0 1]),'.r',MS,10), hold on
plot(y2([0 1]),y2p([0 1]),'.r',MS,10)
arrowplot(y1,y1p,CO,bvpnl,MS,4)
arrowplot(y2,y2p,CO,bvpnl,MS,4), hold off
axis([-5 5 -20 20]), xlabel('$y$',FS,10), ylabel('$y''$',FS,10)
text(1.2,-3,'Fig.~16.2',FS,9)
text(-2.5,9.5,'Fig.~16.3',FS,9)
title(['Fig.~16.6.~~Both solutions to (16.2) ' ...
       'in the phase plane'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Each solution makes a half-circuit clockwise around the origin,
% one on the right starting at $(y,y') = (1,0.7258)$, and
% the other on the left starting at $(y,y') = (1,-13.4074)$.
% By considering this picture we may be tempted to conjecture that in fact,
% (16.2) has infinitely many additional solutions besides these:
% at least two that wind around exactly once,
% two that wind around 1 1/2 times, and so on (Exercise~16.3).
% As the winding
% numbers increase, so do the amplitudes.  For example, here
% is the next solution in the sequence, which we obtain by
% starting with the initial guess $y(x) = 5\sin(2\pi x)$.
% In the phase plane, this corresponds to a large oval winding
% around one full revolution clockwise, beginning at
% $(y,y') = (1,38.8855)$ (not shown).
% </latex>
N.init = 5*sin(2*pi*x);
y3 = N\0; plot(y3,CO,bvpnl); y3p = diff(y3);
hold on, plot(N.init,'--r'), hold off
title(['Fig.~16.7.~~A third solution ' ...
       ' with the corresponding initial guess'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% Whereas nonlinear BVPs can be rather delicate, nonlinear
% IVPs are relatively straightforward.  With this in mind, it is interesting
% to solve (16.2) as an initial value problem, fixing the left-hand
% value at $y(0) = 1$ and taking various choices $y'(0) = a$ in
% the range $[-50,50\kern .3pt]$ for the left-hand derivative value,
% $$ y'' = - y^3, \quad x\in [\kern.3pt 0,1], ~ y(0) = 1, ~ y'(0) = a. 
% \eqno (16.3) $$
%%
% <latex> \vskip -2em </latex>
N.rbc = [];
for a = -50:5:50
  N.lbc = [1;a]; y = N\0; plot(y,LW,0.7,CO,ivpnl), hold on
end
title(['Fig.~16.8.~~Eq. (16.2) as an IVP ' ...
       'with different initial slopes'],FS,11), hold off
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% If we study this image carefully, we can see that there are
% five choices of initial slope $y'(0) = a$ that lead to the condition 
% $y(1)= 1$ being satisfied at the right.  Regarding
% $y(1)$ as a function of $a$, let us plot
% this function over the given range $a\in [-50,50\kern .3pt]$.
% </latex>
fa = @(a) chebop(@(x,y) diff(y,2)+y^3,[0,1],[1;a],[])\0;
f = chebfun(@(a) feval(fa(a),1),[-50,50],120);
plot(f,'k')
title(['Fig.~16.9.~~$y(1)$ as a function of ' ...
       'the initial slope $a = y''(0)$'],FS,11)
xlabel('initial slope $y''(0) = a$',FS,9,IN,LT)
ylabel('$y(1)$',FS,10,IN,LT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Solutions to the BVP (16.2) correspond to points where this curve
% takes the value $1$.
% </latex>
r = roots(f-1); hold on, plot(r,f(r),'.r',MS,12), hold off
title('Fig.~16.10.~~Five solutions of the BVP (16.2)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Five solutions appear in this range,
% </latex>
r'

%%
% <latex>
% \noindent
% and we have seen three of them already,
% </latex>
[y1p(0) y2p(0) y3p(0)]

%%
% <latex>
% \noindent
% In numerical analysis, the technique of varying a 
% slope at one end of an interval so as to satisfy a boundary condition at
% the other end is known as the {\em shooting method.}
% </latex>

%%
% <latex>
% Let us turn to a more celebrated example of nonuniqueness,
% one involving just two solutions rather than an infinite set,
% the {\em Bratu equation:}
% $$ y'' + 3\kern -.7pt  \exp(y) = 0 ,
% \quad x\in [\kern .3pt 0,1], ~ y(0) = y(1) = 0 . \eqno (16.4) $$
% Here we show two solutions on a single plot, one resulting
% from Chebfun's default initial guess (the zero function)
% and the other from the alternative initial guess $8(x-x^2)$.
% </latex>
N = chebop(0,1); N.op = @(x,y) diff(y,2) + 3*exp(y);
N.lbc = 0; N.rbc = 0; y1 = N\0; plot(y1,CO,bvpnl)
hold on, x = chebfun('x',[0 1]); N.init = 8*(x-x.^2);
y2 = N\0; plot(y2,CO,bvpnl), ylim([-.5 2.5]), hold off
title(['Fig.~16.11.~~Two solutions of the ' ...
       'Bratu equation (16.4)'],FS,11)
%%
% <latex> \vskip 1.02em </latex>
%%
% <latex>
% \noindent
% These are the only two solutions to this problem, as we can
% explain (though not quite rigorously prove without some
% more work) by the same method of shooting as before.
% </latex>
N.rbc = [];
for a = -8:2:16
  N.lbc = [0;a]; y = N\0; plot(y,LW,0.7,CO,ivpnl), hold on
end
title(['Fig.~16.12.~~The Bratu eq.\ as an IVP ' ...
       ' with different initial slopes'],FS,11)
ylim([-5 5]), hold off
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% In Chapter 18 we
% shall consider the behavior of equation (16.4) as a parameter
% is varied, namely the coefficient that here takes the value $3$.
% </latex>

%%
% <latex>
% We close this chapter with a final pair of examples
% of nonlinear BVPs with multiple solutions.
% In each case, shooting is a good method to
% explore the variety of solutions, at least when the
% coefficient multiplying the highest derivative
% is not too small.
% Another method of investigation for such problems is
% {\em path-following,} to be considered in Chapter~18.
% Besides these approaches, a further technique is the idea of
% {\em deflation} (Exercise 16.4).
% </latex>

%%
% First we consider the inhomogeneous equation
% $$ \varepsilon y'' + y + y^2 = 1 , \quad x\in [-1,1], ~
% y(\pm 1) = 0. \eqno (16.5) $$
% For $\varepsilon=0.2$
% this problem has four solutions, two symmetric and two asymmetric.
% Higher values of $\varepsilon$ would have just two symmetric
% solutions; the asymmetric ones emerge in a pitchfork
% bifurcation (see next chapter) at a critical value of $\varepsilon$.
N = chebop(-1,1); N.lbc = 0; N.rbc = 0;
N.op = @(x,y) 0.2*diff(y,2) + y + y^2; x = chebfun('x');
N.init = x.^2-1; y1 = N\1; plot(y1,CO,bvpnl)
ylim([-2.6 2.6]), hold on
N.init = 1-x.^2; y2 = N\1; plot(y2,CO,bvpnl)
N.init = sin(pi*x); y3 = N\1; plot(y3,CO,bvpnl)
N.init = -sin(pi*x); y4 = N\1; plot(y4,CO,bvpnl), hold off
title(['Fig.~16.13.~~Four solutions of (16.5) ' ...
       'with $\varepsilon=0.2$'],FS,11,IN,LT)
ylim([-2.5,3])
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% Equation (16.5) is a constant-coefficient variant of an 
% equation known as the {\em Carrier equation,}
% $$ \varepsilon y'' +2(1-x^2)y + y^2 = 1, \quad x\in [-1,1], ~
% y(\pm 1) = 0. \eqno (16.6) $$
% For $\varepsilon=0.2$ this also has four solutions, of similar structure.
% </latex>
N = chebop(-1,1); N.lbc = 0; N.rbc = 0;
N.op = @(x,y) 0.2*diff(y,2) + 2*(1-x^2)*y + y^2;
N.init = x.^2-1; y1 = N\1; plot(y1,CO,bvpnl)
ylim([-2.5 2.5]), hold on
N.init = 1-x.^2; y2 = N\1; plot(y2,CO,bvpnl)
N.init = sin(pi*x); y3 = N\1; plot(y3,CO,bvpnl)
N.init = -sin(pi*x); y4 = N\1; plot(y4,CO,bvpnl), hold off
title(['Fig.~16.14.~~Four solutions of Carrier eq.\ (16.6) ' ...
       'with $\varepsilon=0.2$'],FS,11,IN,LT)
ylim([-2.5 3])
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \begin{center}
% \hrulefill\\[1pt]
% {\sc Application: sending a spacecraft to a destination} \\[-3pt]
% \hrulefill
% \end{center}
% </latex>

%%
% <latex>
% There may be several ways to get a spacecraft from A to B,
% even in a specified time interval.  For example, suppose
% the sun is fixed at the origin in the \hbox{$x$-$y$}~plane and a
% spacecraft starts at position $(-2,1)$.  Starting from the initial velocity
% $(0.7,0.7)$, here is the orbit up to the time $T=6$.  We have
% solved this as a second-order IVP in two variables $x$ and $y$ defining
% the position of the spacecraft.
% </latex>
T = 6; N = chebop(0,T); x0 = -2; y0 = 1; u0 = 0.7; v0 = 0.7;
N.op = @(t,x,y) [diff(x,2) + x/(x^2+y^2)^1.5; ...
                 diff(y,2) + y/(x^2+y^2)^1.5];
N.lbc = @(x,y) [x-x0; y-y0; diff(x)-u0; diff(y)-v0];
orange = [1 .7 0];
[x1,y1] = N\0; plot(0,0,'.',CO,orange,MS,30), hold on, grid on
h1 = arrowplot(x1,y1,CO,ivpnl,LW,1.2,YS,2.5);
axis([-5 5 -1.2 2.8]), set(gca,XT,-4:2:4,YT,0:2:2)
title(['Fig.~16.15.~~A spacecraft trajectory starting ' ...
       'at $(x,y) = (-2,1)$'],FS,11,IN,LT)
%%
% <latex> \vskip 1.02em </latex>

%%
% Of course, actual spacecraft are not simply launched into free
% flight from an initial position and velocity; they have rockets
% and make adjustments along the way.  We are considering
% a simplified problem.

%%
% <latex>
% Now suppose our goal is to choose the initial conditions
% so that the spacecraft will be at 
% position $(2,1)$ at time $T$.
% This is a problem of the kind mentioned in the first footnote 
% of Chapter 5: a BVP in $t$ rather than
% $x$.  If we start from the initial guess
% of a straight line orbit from $(-2,1)$ to $(2,1)$,
% Chebfun finds a solution.
% </latex>
xT = 2; yT = 1;
N.lbc = @(x,y) [x-x0; y-y0]; N.rbc = [xT; yT];
t = chebfun('t',[0 T]); N.init = [-2+2*t/3; 1+0*t];
[x2,y2] = N\0; delete(h1)
arrowplot(x2,y2,CO,bvpnl,LW,1.2,YS,2.5)
title('Fig.~16.16.~~Reaching $(2,1)$ at $T=6$',FS,11,IN,LT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Here are the corresponding initial velocities.
% </latex>
u20 = deriv(x2,0); v20 = deriv(y2,0); disp([u20 v20])

%%
% <latex>
% \noindent
% There is another solution, however, following a longer orbit the
% other way around.  Its average speed will be greater,
% and it will pass closer to the sun.  Here is this second solution plotted
% as a dashed line, obtained by starting from another initial guess.
% </latex>
N.init = [-2+2*t/3; 1-2*sin(pi*t/T)];
[x3,y3] = N\0; arrowplot(x3,y3,LS,'--',CO,bvpnl,LW,1.2,YS,2.5)
title('Fig.~16.17.~~Another solution',FS,11,IN,LT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The vector of initial velocities now points down instead of up.
% </latex>
u30 = deriv(x3,0); v30 = deriv(y3,0); disp([u30 v30])

%%
% <latex>
% \noindent
% Both the solutions we have just found are arcs of ellipses, as the reader
% may confirm by extending the orbits to a later time.  A good
% choice is $T = 27.4$.
% </latex>

%%
% This 2-body problem belongs to
% Newtonian mechanics, and it would have given little difficulty to 
% scholars of earlier centuries.  It is well known that any orbit
% will be an ellipse or a hyperbola, or a parabola in the borderline
% case.  In the present example we got ellipses rather
% than hyperbolas because the number $T$ was sufficiently large.  
% If the spacecraft needs to get from A to B faster, for example in
% $T=2$ time units, then the simplest solution is a nearly-straight
% trajectory that begins pointed almost at $B$, like firing a bullet.
% This time the spacecraft has more than enough kinetic energy to
% escape the sun's gravitational field, and the orbit is an arc of
% a hyperbola rather than an ellipse.
T = 2; N.domain = [0 T];
N.lbc = @(x,y) [x-x0; y-y0]; N.rbc = [xT; yT];
t = chebfun('t',[0 T]); N.init = [-2+2*t; 1+0*t];
hold off, plot(0,0,'.',CO,orange,MS,30), hold on, grid on
[x4,y4] = N\0; arrowplot(x4,y4,CO,bvpnl,LW,1.2,YS,2.5)
axis([-5 5 -1.2 2.8]), set(gca,XT,-5:5,YT,-1:3)
title('Fig.~16.18.~~Faster hyperbolic trajectory: $T=2$',FS,11,IN,LT)

%%
% <latex>
% \noindent
% Again there is a second solution,
% also a hyperbola, that passes the other way around.
% </latex>
N.init = [-2+2*t; 1-1.5*sin(pi*t/T)];
[x5,y5] = N\0; arrowplot(x5,y5,LS,'--',CO,bvpnl,LW,1.2,YS,2.5)
title('Fig.~16.19.~~Another fast trajectory: $T=2$',FS,11,IN,LT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The initial velocities are much greater than before.
% In space travel, a speedy journey may cost a lot of energy!
% </latex>
u40 = deriv(x4,0); v40 = deriv(y4,0); disp([u40 v40])
u50 = deriv(x5,0); v50 = deriv(y5,0); disp([u50 v50])

%%
% As we mentioned, the calculations above are classical.
% When a third body is introduced, however, everything changes
% and only numerical computations are available.  For example,
% suppose that instead of one
% sun at $(0,0)$ we have two fixed stars, one at $(0,0)$ and the other at $(0,1)$.
% All your experience with elliptical orbits is now irrelevant. 
% To see a little of the complexity, let us consider solutions up
% to time $T=12$ beginning with $x'(0)=0.8$ with
% $y'(0) = 0.549$ and $y'(0) = 0.4$.  We see the
% beginnings of two orbits that will remain bounded but
% are certainly not ellipses.
T = 12; N = chebop(0,T); u0 = .8;
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-x0; y-y0; diff(x)-u0; diff(y)-0.549]; hold off
[x6,y6] = N\0; subplot(1,2,1)
plot([0 0],[0 1],'.',CO,orange,MS,20), hold on, grid on
arrowplot(x6,y6,CO,ivpnl,LW,1,YS,.8);
axis([-4 4 -3.5 4.5]), set(gca,XT,-4:2:4,YT,-4:2:4)
N.lbc = @(x,y) [x-x0; y-y0; diff(x)-u0; diff(y)-0.40]; hold off
[x7,y7] = N\0; subplot(1,2,2)
plot([0 0],[0 1],'.',CO,orange,MS,20), hold on, grid on
arrowplot(x7,y7,CO,ivpnl,LW,1,YS,.8,MS,5);
axis([-4 4 -3.5 4.5]), set(gca,XT,-4:2:4,YT,-4:2:4)
title('Fig.~16.20.~~Orbits about two fixed suns',FS,11,HA,RT)
%%
% <latex> \vskip 1.02em </latex>

%%
% Here are the same orbits extended to time $T=140$.  On the left we see
% apparent periodicity.  This is not typical; it results from the particular choice
% $v'(0) = 0.549$.
N.domain = [0 140]; clf
N.lbc = @(x,y) [x-x0; y-y0; diff(x)-u0; diff(y)-0.549]; hold off
[x6long,y6long] = N\0; subplot(1,2,1)
arrowplot(x6long,y6long,CO,ivpnl,LW,.7,YS,.8,MS,5), hold on, grid on
axis([-4 4 -3.5 4.5]), set(gca,XT,-4:2:4,YT,-4:2:4)
plot([0 0],[0 1],'.',CO,orange,MS,20)
N.lbc = @(x,y) [x-x0; y-y0; diff(x)-u0; diff(y)-0.40]; hold off
[x7long,y7long] = N\0; subplot(1,2,2)
arrowplot(x7long,y7long,CO,ivpnl,LW,.7,YS,.8), hold on, grid on
axis([-4 4 -3.5 4.5]), set(gca,XT,-4:2:4,YT,-4:2:4)
plot([0 0],[0 1],'.',CO,orange,MS,20)
title('Fig.~16.21.~~More orbits about two suns',FS,11,HA,RT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Still more complicated orbits can be found by exploring
% values $v'(0)$ such as $0.456$, $0.455$, and $0.454$.
% </latex>

%%
% <latex>
% No solar system, of course, has two stars fixed motionless in space.
% However, our two-star model still has relevance
% to real problems.  First, we could build a system just
% like this involving marbles rolling on a surface with two
% suitably shaped drains in a fixed position (there are demonstrations along these
% lines at some science museums).  Second, there are many binary star systems
% that differ from this only in that the two stars are in orbit around
% each other, which just adds a few more terms
% to the equations.  Third, a single star like our sun may have planets
% orbiting around it, and spacecraft trajectories are often chosen to
% swing close to some of the planets to get to their destination fast
% with less cost in fuel.  The plot of the 2016 movie
% {\em The Martian} turns on such an 
% unexpected choice of orbit.
% </latex>

%%
% <latex>
% \smallskip
% {\sc History.}  We do not know who first focused on the
% phenomenon that a nonlinear BVP can have multiple solutions, but
% certainly one of earliest and most consequential examples of this kind
% concerns the buckling of columns.  Euler published his great
% paper on this subject in 1759 (in French), ``On the strength of
% columns.'' He writes, {\em And so we see that, however small the
% force $F$ acting horizontally, it must always produce a certain
% deflection, which is proportional to $F$ itself.
% But it is not the same when the force
% acts vertically, or if the column must sustain
% a load from above.  At first it seems that such a force, no matter how
% great, could not bend the column: for there is no reason why
% it should bend in one direction rather than another.  But the
% least inequality in the parts of the column, or the
% least stress which it feels from any side, will soon furnish
% a sufficient reason for it to bend in a particular direction.}
% \smallskip
% </latex>

%%
% <latex>
% {\sc Our favorite reference.}  Many features of
% our world, from the fundamental laws of physics to the pattern
% of stripes on a tiger, can be seen as having arisen by the
% process called {\em symmetry breaking}, in which a choice is made
% among a multiplicity of potential solutions.
% For a popular account of such effects see
% Stewart and Golubitsky, {\em Fearful Symmetry: Is God
% a Geometer?,} first published in 1993.
% \smallskip
% </latex>

%%
% <latex>
% \begin{displaymath}
% \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl
% {\sc Summary of Chapter 16.}  Nonlinear ODE BVPs can
% have no solutions, one solution, finitely many solutions,
% or infinitely many solutions.  Some strategies for
% investigating multiple solutions include (a) phase plane
% analysis, (b) shooting, (c) path-following, and 
% (d) deflation.  In general, finding all
% solutions of a nonlinear BVP is difficult.
% \vspace{2pt}}}
% \end{displaymath}
% </latex>

%%
% <latex>\smallskip\small\parskip=1pt\parindent=0pt
% {\em \underline{Exercise $16.1$}.  Fisher equation.}
% The function $y(x)$ satisfies $y''+y-y^2 = 0$
% for $x\in [-1,1]$ with $y(-1) = 1$, $y(1) = 0$.
% {\em (a)} If $y(0)\approx 0.6$, what is $y(0.5)\kern .5pt?$\ \ Plot
% the solution.
% {\em (\kern .7pt b)} If $y(0)\approx -2.5$, what is
% $y(0.5)\kern .5pt?$\ \ Plot the solution.
% {\em (c)} Sketch both of the orbits just described in the phase plane.
% \par
% {\em Exercise $16.2$.  Bounce pass.}
% A ball is thrown from player A to player B, 5 meters away,
% starting and finishing at height 1 meter.  This is an idealized
% ball that travels as a point mass with no air resistance or rotation
% and bounces perfectly with equal angles and speeds of
% impact and rebound.  The pass is a 
% slow one: it takes a full 3 seconds to get from A to B.
% {\em (a)} Assuming the ball does not bounce, sketch its
% trajectory.  You do not need to write any differential equations.
% {\em (\kern .7pt b)} Assuming the ball bounces once, sketch all of its
% possible trajectories.  
% Again you do not need to write any differential equations.
% {\em (c)} Now consider all possible solutions to this BVP,
% with any number of bounces.
% Assume it takes 0.45 seconds for a point mass to fall
% from a height of 1 meter.  Exactly how
% many solutions are there all together?  Sketch them.
% \par
% {\em \underline{Exercise $16.3$}.  Multiple solutions of cubic
% oscillator.}
% Figure 16.6 shows phase plane plots of two of the five solutions of (16.2)
% indicated in Figure 16.10.
% Expand the plot to include all five solutions.  How
% do you think the conjectures stated after Figure 16.6 about
% solutions with various winding numbers should be corrected?
% \par
% {\em \underline{Exercise $16.4$}.  Deflation.}
% {\em (a)} Equation (16.2) can be written $N(y) = 0$, where $N$ is a nonlinear
% operator applying to functions on $[\kern .3pt 0,1]$ with
% boundary values equal to~$1$.  Compute the solution
% to (16.2) shown in Figure 16.3, and
% call this function~$Y$.  Now consider the new nonlinear operator
% $M(y) = N(y)(1 + \|y - Y\|^{-1})$, where $\|\cdot\|$ is the
% 2-norm.  Note that  $M(y) = 0$ for
% $y\ne Y$ if and only if $N(y) = 0$.  Use Chebfun to solve 
% $M(y) = 0$ numerically.  Which
% solution from Figure 16.10 do you get?
% {\em (\kern .7pt b)} This process is automated by
% the Chebfun {\tt deflate} command.  Show that you get the
% same solution with \verb|deflate(N,Y,1,1)|.  Which solution
% do you get with \verb|deflate(N,Y,2,0)|?
% (For information on deflation, see
% Farrell, Birkisson, and Funke, ``Deflation techniques for
% finding distinct solutions of nonlinear partial differential
% equations,'' {\em SIAM Journal on Scientific Computing,} 2015.)
% </latex>
