%% 9. Second-order equations in the phase plane
% <latex>
% \setcounter{page}{96}
% </latex>

%%
% <latex>
% Consider a second-order autonomous ODE, which may be linear or nonlinear,
% $$ y''(t) = f(y,y') . \eqno (9.1) $$
% Given a pair of values $y$ and $y'$ at a particular
% time $t$, this equation tells us the rate of change of $\kern .3pt y'$ at
% $t$, and the rate of change of $y$ itself is by definition~$y'$.
% Thus with (9.1) we know the rates of change of both
% $y$ and $y'$, and we can
% represent this information pictorially by a diagram of the
% {\bf phase plane,} in which
% the horizontal axis represents $y$ and the vertical axis
% represents $y'$.  For example, here is a ``quiver plot''
% of the phase plane
% for the simple harmonic oscillator or linear pendulum
% equation $y'' = -y$.
% </latex>
ODEformats, L = chebop(@(t,y) diff(y,2) + y);
quiver(L,[-2.8 2.8 -1.1 1.1],'k',LW,.3,'scale',.5,'xpts',12,'ypts',8)
title(['Fig.~9.1.~~Phase plane for $y''''=-y$' ...
              ' (simple harmonic motion)'],FS,11,IN,LT)
axis equal, ylim([-1.1 1.1])
hold on, plot(0,0,'.k',MS,8), hold off
xlabel('$y$',FS,10,IN,LT), ylabel('$y''$',FS,11,IN,LT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% At each point $(y,y')$, the arrow shows the direction and
% magnitude of the vector $(y',y'')$.
% Solutions to the ODE will correspond to trajectories
% following the arrows around the plane.
% Note that for any equation (9.1), regardless of $f$, the arrows will
% always point rightward in the upper half-plane and leftward in 
% the lower half-plane.
% The black dot at $y = y'= 0$ marks a 
% {\em fixed point} (or {\em steady state}) of this equation, which
% means, a point at which $y' = y'' = 0$.
% </latex>

%%
% <latex>
% An image like Figure 9.1, while accurate, is sometimes not
% as compelling as one in which all the arrrows are
% set to have the same length, giving a plot of
% a {\em direction field.}  This can be done with the
% {\tt quiver} option \verb|'normalize'|.
% </latex>
quiver(L,[-2.8 2.8 -1.1 1.1],'k',LW,.3,'scale',.4,'normalize',1,'xpts',12,'ypts',8)
title('Fig.~9.2.~~Same but with arrows all of the same length',FS,11,IN,LT)
axis equal, ylim([-1.1 1.1])
hold on, plot(0,0,'.k',MS,8), hold off
xlabel('$y$',FS,10,IN,LT), ylabel('$y''$',FS,10,IN,LT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% There is nothing more to the dynamics of an autonomous
% ODE than its phase plane.  If an
% initial point $(y,y')$ is specified, then all of the future 
% trajectory is determined simply by the vector
% field.\footnote{To be precise, this requires the solutions
% to be unique, which will be assured if $f$ satisfies a
% condition of Lipschitz continuity.  See Chapter 11.}
% This is the power of phase plane analysis: it reduces dynamics
% to geometry.
% </latex>

%%
% <latex>
% For example, let us now specify a\/ $\kern .5pt t\kern .5pt$ domain and an initial condition
% for the linear pendulum,
% $$ y''= -y, \quad t\in [\kern .5pt 0, 1.8\pi],~
% y(0) = 0, ~ y'(0)=1, \eqno (9.2) $$
% or in Chebfun,
% </latex>
L.domain = [0,1.8*pi]; L.lbc = [0;1];
%%
% <latex>
% \noindent
% Here we superimpose the trajectory corresponding to the
% solution of (9.2) on the vector field just displayed.
% The trajectory begins at the top of the unit
% circle with initial velocity~$1$
% and initial acceleration~$0$, so the curve is
% oriented horizontally to the right.  As $y$ increases, $y''$ becomes 
% negative, and the curve begins to bend around, describing
% a clockwise circle.
% Since the time interval runs $90\%$ of the way to $2\pi$,
% the circle is $90\%$ complete.
% </latex>

y = L\0; hold on, arrowplot(y,diff(y),CO,ivp,YS,3.5)
title(['Fig.~9.3.~~A trajectory added to the plot, ' ...
       'the solution of (9.2)'],FS,11), hold off
%%
% <latex> \vskip 1.02em </latex>

%%
% Let us dispense with the quiver arrows and collect
% trajectories with different initial values on a single
% plot.  The
% next picture takes $b=0.5,0.8,\dots, 2$, producing six circles
% of corresponding radii.
for b = .5:.3:2
  L.lbc = [0;b]; y = L\0;
  arrowplot(y,diff(y),CO,ivp,MS,4,YS,3), hold on
end
axis equal, ylim([-2.2 2.2])
plot(0,0,'.k',MS,8)
title('Fig.~9.4.~~Six trajectories for (9.2)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% The phase plane can be informative about boundary-value
% problems, too.  For example, suppose we change (9.2) to
% $$ y''= -y, \quad t\in [\kern .5pt 0, 3],
% ~ y(0)=1, ~ y(4) = -1.5. \eqno (9.3) $$
% Here is the solution, shown as a seventh trajectory added
% to the previous plot (without an arrowhead since this
% is a BVP).
L = chebop(0,4); L.op = @(t,y) diff(y,2) + y;
L.lbc = 1; L.rbc = -1.5; y = L\0; plot(y,diff(y),CO,bvp)
hold off, title('Fig.~9.5.~~Solution of BVP (9.3)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% Equations (9.2) and (9.3) describe a linear pendulum without 
% damping.  Following Chapters 4 and 8, let us now add
% some damping in the form of a small multiple of $y'$,
% $$ y''= -y-\varepsilon \kern .3pt y', \quad t\in [\kern .5pt 0, 1.8\pi],
% ~ y(0)=0, ~ y'(0) = b, \eqno (9.4) $$
% with $\varepsilon=0.2$.
% If we follow the 
% same six trajectories as before, we see that they now lose amplitude
% as they evolve, moving towards the fixed point
% $y = y'= 0$ as $t\to\infty$.
% This is a {\em stable} fixed point since all nearby trajectories
% stay close to it; we shall consider such definitions
% systematically in Chapter~15.\ \ The
% plot also shows the solution to the BVP variant of
% (9.4) with the same boundary conditions $y(0)=1$,
% $y(4) = -1.5$ as in (9.3).
% Note that the solution starts from a greater velocity $y'$
% than before, necessary to compensate for
% the damping while still matching the boundary conditions.
% </latex>
L = chebop(0,1.8*pi); L.op = @(t,y) diff(y,2) + 0.2*diff(y) + y;
for b = .5:.3:2
  L.lbc = [0;b]; y = L\0;
  arrowplot(y,diff(y),CO,ivp,MS,3,YS,3), hold on
end
axis equal, ylim([-2.2 2.2])
L = chebop(0,4); L.op = @(t,y) diff(y,2) + 0.2*diff(y) + y;
L.lbc = 1; L.rbc = -1.5; y = L\0; plot(y,diff(y),CO,bvp)
plot(0,0,'.k',MS,8), hold off
title('Fig.~9.6.~~Damped linear pendulum (9.4)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% We stopped the flow at $t=1.8\pi$, but of course this is not necessary.
% Here is a single trajectory carried to $9\pi$.
% This represents a linear pendulum swinging back and forth $4~1/2$ times,
% losing amplitude as it swings.
% </latex>
L = chebop(0,9*pi); L.op = @(t,y) diff(y,2)+0.2*diff(y)+y;
L.lbc = [0;2]; y = L\0;
plot(y,diff(y),CO,ivp), axis equal, ylim([-2.2 2.2])
hold on, plot(0,0,'.k',MS,8), hold off
title(['Fig.~9.7.~~A trajectory for (9.4) ' ...
       'carried to longer time'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% We mentioned that in
% the phase plane, the position of a trajectory at any time
% $t$ determines its future, which consists of following the unique path
% through that point.  This uniqueness property
% depends on the equation being autonomous.  If
% we solve a nonautonomous problem, such as
% $$ y''= -y-0.2\kern .3pt y'-2\cos(2\kern .3pt t), \quad t\in [\kern .5pt 0, 1.8\pi],
% ~ y(0)=0, ~ y'(0) = b, \eqno (9.5) $$
% we can plot the solution in the $y$-$y'$ plane, like this (now
% with different colors to help distinguish the curves),
% </latex>
L = chebop(0,1.8*pi);
L.op = @(t,y) diff(y,2)+0.2*diff(y)+y+2*cos(2*t);
for b = .5:.3:2
  L.lbc = [0;b]; y = L\0; arrowplot(y,diff(y),LW,.5,MS,3), hold on
end
axis equal, ylim([-3 2.4]), plot(4.5,0,'Xr',MS,30,LW,2)
title(['Fig.~9.8.~~Inhomogeneous problem --- not truly '...
        'a phase plane!'],FS,11), hold off
%%
% <latex> \vskip 1.02em </latex>
%%
% <latex>
% \noindent
% However, the point of planar analysis has
% been lost since the future of a trajectory is not
% determined by its current position in the plane, and the
% curves cross each other.  The plot of trajectories
% has become merely a plot of trajectories, no longer
% an encapsulation of the dynamics of the system.
% Such plots can be interesting and informative; see for example
% Exercise 9.2 and the Duffing oscillator in Appendix~B.\ \ They
% usually don't belong to phase plane analysis, however.
% An alternative in such cases is to include the time variable in
% the plot on an additional axis.  Here, for example, we repeat the
% calculation on the longer interval $[\kern .3pt0 ,40\kern .3pt]$,
% for a single trajectory, and plot the result in $t$-$y$-$y'$ space.
% </latex>
L.domain = [0,40]; y = L\0;
t = chebfun('t',[0 40]); plot3(t,y,diff(y),CO,ivp,LW,.5), view(-48,27)
axis([0 40 -2.5 2.5 -3 2])
axis square, xlabel('$t$',FS,10), ylabel('$y$',FS,10), zlabel('$y''$',FS,10)
title('Fig.~9.9.~~Trajectory in $t$-$y$-$y''$ space',FS,11)
set(gcf,'position',[380 220 540 300])
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% Returning to (9.4), in the language of Chapter 4, the damping
% with $\varepsilon = 0.2$ is {\em subcritical\/}: from (4.11) we see that
% the critical damping parameter for this equation is $\varepsilon=2$.
% Here are images of subcritical, critical,
% and supercritical trajectories with
% $\varepsilon = 1,2,4$, now for $t\in [\kern .3pt 0, 1.5\pi]$.
% Note that it is the
% critically damped trio of trajectories that finishes closest
% to the origin, confirming that $\varepsilon=2$ gives 
% the most effective damping.
% </latex>
close all, L.domain = [0,1.5*pi];
subplot(1,3,1), L.op = @(t,y) diff(y,2)+diff(y)+y;
for b = 1:-.2:.6
  L.lbc = [0;b]; y = L\0;
  arrowplot(y,diff(y),MS,3,YS,.8), hold on
end
axis equal, ylim([-0.4 1.25])
pp = get(gca,'pos'); pp(1) = pp(1)+.03; set(gca,'pos',pp)
text(.21,1.09,'subcritical, $\varepsilon=1$',FS,9,HA,CT)
hold on, plot(0,0,'.k',MS,8), hold off
subplot(1,3,2), L.op = @(t,y) diff(y,2)+2*diff(y)+y;
for b = 1:-.2:.6
  L.lbc = [0;b]; y = L\0; arrowplot(y,diff(y),MS,3,YS,.8), hold on
end
axis equal, ylim([-0.4 1.25])
pp = get(gca,'pos'); pp(1) = pp(1)-.0003; set(gca,'pos',pp)
set(gca,YT,0:.5:1,YTL,{' ',' ',' '})
text(.21,1.09,'critical, $\varepsilon=2$',FS,9,HA,CT)
hold on, plot(0,0,'.k',MS,8), hold off
title(['Fig.~9.10.~~Subcritical, critical, and ' ...
       'supercritical damping'],FS,11,HA,CT)
subplot(1,3,3), L.op = @(t,y) diff(y,2)+4*diff(y)+y;
for b = 1:-.2:.6
  L.lbc = [0;b]; y = L\0; arrowplot(y,diff(y),MS,3,YS,.8), hold on
end
axis equal, ylim([-0.4 1.25])
set(gca,YT,0:.5:1,YTL,{' ',' ',' '})
pp = get(gca,'pos'); pp(1) = pp(1)-.03; set(gca,'pos',pp)
text(.12,1.09,'supercritical, $\varepsilon=4$',FS,9,HA,CT)
hold on, plot(0,0,'.k',MS,8), hold off
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% Phase plane analysis becomes particularly interesting
% for nonlinear ODE\kern .5pt s.  For example, the
% first nonlinear equation of this book was the
% van der Pol equation.  With slightly
% different coefficients from (1.2), let us write this as
% $$ y'' + y - \mu (1-y^2)y' = 0, \quad t\in[\kern .5pt 0, 10\kern .5pt],
% ~ y(0) = a, y'(0) = 0. \eqno (9.6) $$
% Here are phase plane plots for $\mu = 0.125$ and $\mu = 1.5$.
% With the weak damping parameter $\mu=0.0125$, the system is
% not far from the linear pendulum, with trajectories winding
% slowly in or out to an asymptotic curve known
% as a {\bf limit cycle}.  With the stronger damping parameter
% $\mu = 1.5$, trajectories converge to the limit cycle much faster,
% and its shape is far from circular, just as the
% van der Pol orbit plotted in Chapter 1 is far from a sine wave.
% </latex>
N = chebop(0,15); clf
for j = 1:2
  subplot(1,2,j), mu = 0.125; if j==2, mu = 1.5; end
  N.op = @(t,y) diff(y,2)-mu*(1-y^2)*diff(y)+y;
  for a = [1 3]
    N.lbc = [a;0]; y = N\0; arrowplot(y,diff(y),LW,.6,MS,4), hold on
  end
  axis equal, ylim([-4.2 4.2]), plot(0,0,'.k',MS,8)
  text(4.5,3.1,['$\mu = ' num2str(mu) '$'],FS,9,HA,RT)
end
subplot(1,2,1)
title('\kern 2.4 in Fig.~9.11.~~Van der Pol equation',FS,11)
%%
% <latex> \vskip 1.02em </latex>

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

%%
% <latex>
% The linear pendulum equation $y'' = -y$ describes oscillation in
% the context of Hooke's Law, where the restoring force is
% proportional to the displacement --- simple harmonic motion.
% When the amplitudes of motion are small,
% this is the right model for a spring or a
% pendulum and for many other vibrating systems, but when
% the amplitudes get bigger, the physics always becomes nonlinear.
% Different problems have different nonlinearities, but
% there is no doubt as to the archetypical problem of this
% kind: it is the {\bf nonlinear pendulum},
% corresponding to an idealized point mass moving in a
% circle at the end of a rigid weightless bar.
% The equation is
% $$ y''= -\sin(y), \eqno (9.7) $$
% where $y(t)$ represents the angle from the vertical in radians
% at time $t$ and constants are set to 1.
% For an entire book on the subject,
% see Baker and Blackburn, {\em The Pendulum,}
% Oxford University Press, 2005.
% </latex>

%% 
% <latex>
% This is a perfect example for phase plane analysis.
% First we draw a quiver plot together with stable and
% unstable fixed points (black and red, respectively).
% There are stable fixed points at $(y,y') = (2\pi j,0)$
% for each integer~$j$ and unstable fixed
% points at the in-between locations $(2\pi(\kern .7pt j+1/2),0)$.
% </latex>
N = chebop(0,1.8*pi);
N.op = @(t,y) diff(y,2)+sin(y); clf
quiver(N,[-3 23 -6 6],'k',LW,.3,'scale',.5,'xpts',12,'ypts',8)
axis equal, xlim([-3 23]), hold on
plot(pi*[0 2 4 6],[0 0 0 0],'.k',MS,10)
plot(pi*[1 3 5 7],[0 0 0 0],'.r',MS,10), hold off
title('Fig.~9.12.~~Nonlinear pendulum (9.7)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% To see more, let us plot some trajectories.
% We start from $y=0$ and
% $y'(0) = b = 1,1.2,\dots, 4$, calculating trajectories
% over the interval $t\in [\kern .5pt 0, 1.8\pi]$.\footnote{To us this
% looks like the hair of Botticelli's Venus; or is it
% the Starbucks logo?}
% </latex>
plot(pi*[0 2 4 6],[0 0 0 0],'.k',MS,10), hold on
plot(pi*[1 3 5 7],[0 0 0 0],'.r',MS,10)
for b = 1:.2:4
  N.lbc = [0;b]; y = N\0;
  arrowplot(y,diff(y),LW,.6,CO,ivp,YS,3,MS,3)
end
axis equal, xlim([-3 23]), hold off
title('Fig.~9.13.~~Nonlinear pendulum (9.7)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Physically, these curves can be interpreted as follows.
% If $y'(0)=b<2$, the pendulum does not have enough
% energy to swing over the top, and the trajectory goes
% around and around forever on a single loop in the
% phase plane, corresponding to the pendulum swinging
% back and forth, with the angle $y(t)$ remaining
% bounded.  The loops are not circles, but they
% approach circles for small amplitude, where the
% distinction between $y$ as in (9.2) and $\sin(y)$ 
% as in (9.7) fades away.
% If $y'(0)=b>2$, on the other hand,
% the pendulum has enough energy to 
% swing over, and $y(t)$ keeps increasing monotonically.
% In the absence of damping, a pendulum that swings over
% once swings over infinitely many times as $t\to\infty$.
% The phase plane is $2\pi$-periodic.
% </latex>

%%
% Notice the trajectory starting at $y'(0) = 2$.  This one appears
% to stop at the $y$ axis, and that is exactly what it does.
% It has just enough energy to fly up toward the vertical 
% configuration, with $y(t)$ approaching $\pi$ as $t\to\infty$, but
% it never quite hits the top for any finite value of $t$.

%%
% Of course a real pendulum is sure to have some losses.
% Here is the same image corresponding to an equation with subcritical
% damping,
% $$ y''= -\sin(y) - \mu y', \quad t\in [\kern .5pt 0, 1.8\pi],
% ~ y(0)=0, ~ y'(0) = b, \eqno (9.8) $$
% with $\mu = 0.1$.
% A curve is also added to the plot corresponding to the
% BVP defined by $y(0)= 0$, $y(1.95\pi) = 20$.
% We see that to reach $y=20$ by the end of
% the time interval, a greater initial speed is needed.
N.op = @(t,y) diff(y,2)+sin(y)+.1*diff(y);
for b = 1:.2:4
  N.lbc = [0;b]; y = N\0;
  arrowplot(y,diff(y),LW,.6,CO,ivp,YS,2,MS,3), hold on
end
plot(pi*[0 2 4 6],[0 0 0 0],'.k',MS,10)
plot(pi*[1 3 5 7],[0 0 0 0],'.r',MS,10)
axis equal, xlim([-3 23])
title('Fig.~9.14.~~Damped nonlinear pendulum (9.8)',FS,11)
N.lbc = 0; N.rbc = 20;
y = N\0; plot(y,diff(y),CO,bvpnl), hold off
%%
% <latex> \vskip 1.02em </latex>
%%
% <latex>
% \noindent
% Now, a trajectory
% with sufficient initial energy swings over, but it
% keeps losing energy as it goes, so it may or may
% not swing over a second time.  We can see more
% if we increase the time interval to $[\kern .5pt 0,8\pi]$.
% The trajectory with initial velocity $y'(0) = 4$ swings
% over four times before eventually winding down to rest.
% As before, the function $f(y,y')$ defining the phase plane
% is $2\pi$-periodic, though none of the individual trajectories
% are periodic.
% </latex>
N = chebop(0,8*pi); N.op = @(t,y) diff(y,2)+sin(y)+.1*diff(y);
plot(pi*[0 2 4 6 8],[0 0 0 0 0],'.k',MS,10), hold on
plot(pi*[1 3 5 7 9],[0 0 0 0 0],'.r',MS,10)
for b = 1:.2:4
  N.lbc = [0;b]; y = N\0; arrowplot(y,diff(y),LW,.6,CO,ivp,MS,3)
end
axis equal, xlim([-3 29]), hold off
title('Fig.~9.15.~~Longer time interval for (9.8)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \medskip
% {\sc History.}  Four hundred years ago, Galileo understood
% the essentials of the linear
% pendulum.  In his {\em Dialogues Concerning Two New
% Sciences\/} the character Salviati
% says, ``As to the times of vibration of bodies suspended by
% threads of different lengths, they bear to each other the
% same proportion as the square roots of the lengths of
% the threads; or one might say the lengths are to each other
% as the squares of the times; so that if one wishes to make
% the vibration-time of one pendulum twice that of another,
% he must make its suspension four times as long.''
% </latex>

%%
% <latex>
% \smallskip
% {\sc Our favorite reference.}  When it comes to phase plane
% analysis, one of the oldest books is also one of the
% nicest: H. T. Davis, {\em Introduction to Nonlinear
% Differential and Integral Equations,} Dover, 1962 (first published
% in 1960).  For example, we enjoy Chapter 10 on ``The phase plane
% and its phenomena,'' Chapter 11 on ``Nonlinear
% mechanics,'' and Chapter 12 on ``Some particular equations.''
% \smallskip
% </latex>

%%
% <latex>
% \begin{displaymath}
% \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl
% {\sc Summary of Chapter 9.}   The phase plane for a
% second-order autonomous homogeneous ODE\/ $y'' = f(y,y')$ is
% the plane with coordinates\/ $y$ and $y'$.  The position of
% a solution at a particular time\/ $t$ determines its future trajectory
% as the unique curve passing through this point.  Plotting
% phase plane trajectories gives a quick way to interpret behavior of
% ODE\kern .5pt s, including nonlinear ones.
% \vspace{2pt}}}
% \end{displaymath}
% </latex>

%%
% <latex>
% \small\parindent=0pt\parskip=1pt
% {\em Exercise $9.1$.  $y'' = y$.}
% Draw phase plane trajectories or
% quiver arrows for the equation $y'' = y$ (either
% by hand or with the computer).
% Where in this image does the solution of Fig.~5.3 lie?
% (It comes very close to following the {\em stable manifold} of
% the fixed point for 30 time units and then the
% {\em unstable manifold} for a further 30 times units; see 
% Chapter~15.)
% \par
% {\em \underline{Exercise $9.2$}.  Moving slowly in
% the phase plane.}
% Consider the ODE $y'' + y = y^2$.
% {\em (a)} What are the fixed points?
% {\em (\kern .7pt b)} Draw a quiver plot, choosing axes to make
% the plot as informative as possible.  Based on this
% information, draw a sketch by hand of the key points
% of the phase plane dynamics.
% Describe qualitatively what orbits will remain bounded,
% and how they will behave; also what orbits will diverge
% to $\infty$, and how they will behave.
% {\em (c)} Find two distinct solutions that satisfy
% $y(0) = y(10) = 2$.  Plot them both as functions of $t$ and
% in the phase plane.
% What are the values of $y'(0)$ for these two solutions?
% \par
% {\em \underline{Exercise $9.3$}.  Region-filling orbits.}
% (Adapted from Davis, {\em Introduction to Nonlinear
% Differential and Integral Equations,} \S 10.4.)
% {\em (a)} Use the method of undetermined coefficients to find the
% analytical solution of the IVP $y''+2\kern .3pt y =
% -2\cos(2\kern .3pt t)$, $y(0)=1$, $y'(0)=\sqrt{2}$. 
% {\em (\kern .7pt b)} This solution has a curious property in the phase plane: the
% curve eventually gets arbitrarily close to every point in a
% near-elliptical region. Plot the solution over $t\in[\kern .3pt 0,100\kern .3pt]$ to see the
% effect. (This is fun to watch using {\tt comet} as well.) 
% {\em (c)} The solution just considered is nonperiodic.  Periodic
% solutions are obtained from the same initial conditions, on the
% other hand, if $\cos(2\kern .3pt t)$ is replaced
% by $\cos(k\kern .3pt t)$ for certain
% values of $k>0$.  Find the smallest three such values analytically,
% and produce the corresponding plot for $t\in[\kern .3pt 0,100\kern .3pt]$.
% \par
% {\em \underline{Exercise $9.4$}.  Period of the nonlinear pendulum.} 
% The code
% \begin{verbatim}
%   y = @(s) chebop(@(y) diff(y,2)+sin(y),[0 15-5*log(pi-s)],[s;0],[])\0;
%   T = @(s) 2*min(diff(roots(y(s))));
% \end{verbatim}
% produces two anonymous functions $y$ and $T$ applicable for
% values $s\in (0,\pi)$.
% Explain what these functions compute and how they do it.
% Make a plot of $T(s)$ for $s\in[\kern .3pt 0.1,3.14]$
% and explain its principal features.
% Explain in a general way (not necessarily with 
% mathematical details) the role of the quantity $15-5\log(\pi-s)$.
% </latex>
