%% 21. Into the complex plane

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

%%
% Most of the familiar mathematical functions are defined for
% complex arguments as well as real ones.  For example, here are the 
% Taylor series of two functions about the point $t_0=0$:
% $$ e^{\kern .3pt t} = 1 + t + {t^2\over 2} + {t^3\over 3!} + \cdots, \quad 
% {1\over 1+t^2} = 1 - t^2 + t^4 - t^6 + \cdots. $$
% To define the value of a function when $t$ is complex,
% we can use the same series, so long as $t$ is close
% enough to $t_0$ for convergence.  

%%
% <latex>
% We say that a function $y(t)$ is {\bf analytic} at a point
% $t_0$ if it has a Taylor series at~$t_0$ that converges to $y$ in some
% neighborhood of $t_0$.  If $y$ is initially defined just for real values,
% then the neighborhood will be an interval $(a,b)$ with
% $a < t_0 < b$.  Once we've got a convergent power series, it can
% be used to define values in the complex plane too,
% and this is the process called
% {\bf analytic continuation}.  It is a basic fact
% of complex variables that every power series has a
% {\bf radius of convergence} $r\in [\kern .3pt 0,\infty]$ with the property
% that the series converges for all real or complex
% numbers $t$ with $|\kern .5pt t-t_0|<r$ and
% diverges for all $t$ with $|\kern .5pt t-t_0|>r$.  Thus a power
% series always converges inside a disk, the
% {\bf disk of convergence,} which is the largest disk
% centered at $t_0$ inside which the function is analytic.  
% </latex>

%%
% <latex>
% For $y(t) = e^{\kern .3pt t}$, the disk of convergence has radius
% $r=\infty$, regardless
% of the choice of $t_0$, because the exponential function is analytic
% for all values of $t$.  A function $y(t)$ that 
% is analytic for all real and complex values
% of $t$ is said to be {\bf entire}.  For $1/(1+t^2)$ and
% $t_0=0$, we have $r=1$, because the function has 
% singularities at $t=\pm i$ (poles), where the value blows
% up to $\infty$.  For $t_0=1$ with the same function, we
% would have $r=\sqrt 2$.
% </latex>

%%
% <latex>
% It is a general idea of mathematics that one can learn
% things about a function by examining its singularities in the
% complex plane.  For example, how would you explain to 
% a student in a calculus class {\em why\/} the function
% $y(t) = 1/(1+t^2)$ has a Taylor series 
% that converges just for $t\in (-1,1)$?  The answer
% is the presence of those singularities at $\pm i$, and
% if your student doesn't know about the complex plane, it
% is not clear how you could really give a satisfactory explanation.
% </latex>

%%
% <latex>
% Here is another example.  The function $y(t) = \tanh(10\kern .7pt  t)$
% makes a rapid transition from $-1$ to $+1$ near $t=0$.
% </latex>
ODEformats, y = chebfun('tanh(10*t)'); plot(y), ylim([-1.5 1.5])
title('Fig.~21.1.~~A function $y(t)$ with a rapid transition',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% How does this transition come about?
% One way to understand it is to note that~$y$ has poles 
% in the complex plane quite near to $t=0$, at $\pm \pi i/20,$ $\pm 3\pi i/20,$ 
% $\pm 5\pi i/20,\dots.$ 
% As $t$ increases along the real axis, $y$
% changes quickly as it passes these special points.
% A contour plot of $|y(t)|$ for $t$ ranging over a rectangular region
% in the complex plane reveals the four poles of $y$
% closest to the origin. 
% </latex>
x = linspace(-1.2,1.2,201); y = linspace(-0.5,0.5,101);
[xx,yy] = meshgrid(x,y); zz = xx + 1i*yy;
ff = tanh(10*zz); levels = [1 2 4];
contour(x,y,abs(ff),levels,'k',LW,.7), axis equal
hold on, plot(pi*.05i*(-3:2:3),'.r',MS,8), hold off
title(['Fig.~21.2.~~In the complex $t$-plane, ' ...
       'there are poles nearby'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% If the poles were ten times closer, $y$ would make its transition ten times faster.
% </latex>

%%
% For another example, what can we make of the function
% $$ y(t) = -\log (1.1 + \cos(\pi t))?  \eqno (21.1) $$
% A plot reveals spikes near odd integer values of $t$.
f = chebfun('-log(1.1+cos(pi*t))',[-4 4]); plot(f), ylim([-3 3])
title(['Fig.~21.3.~~Another function $y(t)$ ' ...
       'with regions of rapid transition'],FS,11)
%%
% <latex> \vskip 1.02em </latex>
%%
% <latex>
% \noindent
% As before, these points of rapid change are
% associated with nearby singularities in the complex plane, which
% lie in pairs just above and below each odd integer.
% </latex>
x = linspace(-4,4,201); y = linspace(-1.5,1.5,101);
[xx,yy] = meshgrid(x,y); zz = xx + 1i*yy;
ff = -log(1.1+cos(pi*zz)); levels = 2.5:.5:7;
contour(x,y,abs(ff),levels,'k',LW,.7), axis equal
hold on, plot(acos(1.1)/pi+(-3:2:3),'.r',MS,8)
plot(-acos(1.1)/pi+(-3:2:3),'.r',MS,8), hold off
title(['Fig.~21.4.~~Now the nearby complex ' ...
       ' singularities are branch points'],FS,11)
%%
% <latex> \vskip 1.02em </latex>
%%
% <latex>
% \noindent
% This contour plot requires more interpretation than the previous
% one.  Before, each singularity $t_s$ was a 
% {\bf pole,} with $y\sim C/(t-t_s)$ for some constant $C$.
% In this second
% contour plot, however, the two rows of singularities
% are not poles but {\bf branch points}.  A function with a branch
% point is not analytic and single-valued in any punctured
% neighborhood of the point, because if you continue it
% analytically all the way around, you reach a different value from
% the one you started with.
% To make a function $y$ with branch points single-valued and analytic,
% it is customary to restrict the domain by introducing {\bf branch cuts} connecting the
% branch points, along which $y$ is not defined.  Given a function
% with certain branch points, there is no unique
% set of branch cuts that the mathematics forces you to choose.
% Computer programming languages, however, generally 
% take $\log(t)$ and $t^{1/2}$ to have
% branch lines extending along the negative real axis from $0$ to
% $-\infty$, and this determines where the branch cuts lie for
% more complicated functions like (21.1).
% In the figure, the branch lines have ended up
% as vertical rays extending upward to $\infty$ from the
% branch points in the upper half-plane
% and downward to $\infty$ from the branch points in the lower
% half-plane.
% </latex>

%%
% <latex>
% Everything we have said in the last three
% pages about functions in general applies to
% the particular case of solutions of ODE\kern .5pt s.  Many of them make sense
% in the complex plane, and their behavior in the plane can
% both be interesting in its own right and also shed light on
% behavior, such as rapid transitions, on the real axis.
% </latex>

%%
% As the simplest possible example, consider the IVP
% $$ y' = y, \quad y(0) = 1 . \eqno (21.2) $$
% The unique solution is $e^{\kern .3pt t}$, an entire function in the complex
% $t$-plane, and thus for example at $t=i$ we have
% $$ y(i) = \exp(i) = \cos(1) + i\kern .3pt \sin(1)
% \approx 0.5403 + 0.8415i . $$
% We can interpret (21.2) in the complex plane
% in two ways.  One is to
% regard it as a differential equation that applies not
% just for real values of $\kern .3pt t$ but also 
% complex ones.  The value $y(i)$, for example, might
% be determined by integrating the ODE along the imaginary line segment
% extending from $t=0$ to $t=i$.  
% Chebfun doesn't work with complex intervals directly, but
% we can achieve the necessary effect by parametrizing the
% interval $[\kern .3pt 0, i\kern .3pt ]$ as $i s$ for
% $s\in [\kern .3pt 0,1]$.  The equation becomes
% $$ {dy\over ds} = {dt\over ds} {dy \over dt} = i\kern .3pt y, 
% \eqno (21.3) $$
% or in Chebfun, 
L = chebop(0,1); L.lbc = 1;
L.op = @(s,y) diff(y) - 1i*y; y = L\0; y(1)
%%
% <latex> \vskip 1.02em </latex>

%%
% The other interpretation of (21.2)  for complex $t$ is that 
% we may start from the solution $y(t)$ for real $t$, and then
% analytically continue it into the complex plane.  This idea defines the
% same function as before, provided one uses the same branch
% cuts for the analytic continuation as for the ODE solution.

%%
% For example, here we solve (21.2) in the usual manner
% on the real interval
% $[\kern .3pt 0,1]$, and then evaluate the result
% at the complex point $t = i$.
L = chebop(0,1); L.op = @(t,y) diff(y) - y; 
L.lbc = 1; y = L\0; y(1i)
%%
% <latex> \vskip 1.02em </latex>
%%
% <latex>
% \noindent
% Since Chebfun's solutions of ODE\kern .5pt s are numerical rather
% than symbolic, it is entirely appropriate to wonder
% why this worked!
% The reason is that Chebfun's solution function $y$
% is actually a polynomial representation of $\exp(t)$
% on $[\kern .3pt 0,1]$, which also approximates the
% function in a larger region of the complex plane.
% So Chebfun has performed a numerical version of 
% analytic continuation for us.
% As a rule, this trick works reasonably well for complex
% values of $\kern .3pt t$ close to the interval of definition,
% but not for values further out, and certainly not for
% values of $\kern .3pt t$ that lie further from the interval
% than some singularities of $y$.
% </latex>

%%
% Let us illustrate this by considering an example with a singularity.
% The nonlinear IVP
% $$ y' = -2\kern .3pt t\kern .3pt y^2, \quad t\in [-4,4],
% ~~ y(-4) = {1\over 17}  \eqno (21.4) $$
% has the unique solution we examined at the start of the chapter,
% $$ y(t) = {1\over 1+t^2}, $$
% with poles at $t = \pm i$.  The solution looks as it should,
N = chebop(-4,4); N.lbc = 1/17; N.op = @(t,y) diff(y) + 2*t*y^2;
y = N\0; plot(y,CO,ivpnl), ylim([0 1.2])
title(['Fig.~21.5.~~Solution of an ODE ' ...
       'passing near a complex singularity'],FS,11)
%%
% <latex> \vskip 1.02em </latex>
%%
% <latex>
% \noindent
% At $t=i/2$, the correct value of $y$ is $4/3$,
% as Chebfun confirms,
% </latex>
y(1i/2)
%%
% <latex>
% \noindent
% At $t=3i/2$, however, on the far side of the pole at $i$,
% the correct value of $y$ is $-4/5$, whereas
% Chebfun gets a huge incorrect value.
% </latex>
y(3i/2)
%%
% <latex>
% \noindent
% One can never expect a simple evaluation like this to work,
% if the point of evaluation is further from the interval of
% computation than some of the singularities of the
% function.\footnote{This statement can be made mathematically
% precise with the use of {\em Bernstein ellipses,} ellipses
% in the complex $t$-plane whose foci are at the two endpoints
% of the interval of definition of the chebfun.  A chebfun
% will normally give a degree of good approximation in the largest
% region of analyticity bounded by such an ellipse, 
% as discussed in Chapter~8 of
% Trefethen, {\em Approximation Theory and Approximation Practice,}
% SIAM, 2013.}
% </latex>

%%
% <latex>
% There is an alternative approach in Chebfun, however, that can often
% be used for numerical analytic continuation even beyond singularities,
% especially when the singularities are poles.  The Chebfun command
% {\tt aaa} approximates a function not by the usual
% Chebfun method involving polynomials, but by a rational function
% $r(t) = p(t)/q(t)$ of some adaptively determined type $(n,n)$, which means 
% that the numerator and denominator degrees are $\le n$.\footnote{See
% Nakatsukasa, S\`ete, and Trefethen, ``The AAA algorithm
% for rational approximation,'' submitted manuscript, 2016.} For example, the command
% </latex>
[r,pol] = aaa(y,'tol',1e-8);

%%
% <latex>
% \noindent
% constructs a rational approximation to $y$ and evaluates its poles, which 
% match those of $y$,
% </latex>
pol

%%
% <latex>
% \noindent
% This time, the value at $t = 3i/2$ comes out correct.
% </latex>
r(3i/2)

%%
% <latex>
% Although rational approximations will usually not give
% precise information, especially
% since singularities of solutions to ODE\kern .5pt s in the complex
% plane are usually more complicated than just poles, they
% are often good at giving a rough idea of the nature of such
% singularities near the domain of approximation.  For example, in Chapter 13
% we plotted one component of a solution to the Lorenz equations (13.1),
% </latex>
N = chebop(0,5);
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,CO,ivpnl);
title('Fig.~21.6.~~A trajectory of the Lorenz equations',FS,11)
%%
% <latex> \vskip 1.02em </latex>
%%
% <latex>
% \noindent
% You might guess that these spikes are associated with nearby
% singularities in the complex plane, and here are the poles of the
% {\tt aaa} rational approximation.
% Note that the sharpest spikes correspond to
% singularities closest to the real axis.
% </latex>
[r,pol] = aaa(u,'tol',1e-8); 
plot(pol,'.r',MS,10), xlim([0 5]), axis equal
title(['Fig.~21.7.~~Nearby complex singularities ' ...
       'found by rational approximation'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex> 
% One should not imagine that the function plotted in Fig.~21.6 truly
% has poles at all the points shown in Fig.~21.7.  Rather, it is likely that
% it has a singularity --- probably not just a pole --- near the
% point in each group that lies closest to the real axis.
% The rest of the dots are probably lining up along branch cuts.
% </latex>

%%
% <latex>
% Solutions to linear ODE\kern .5pt s have no singularities unless the leading-order
% coefficient passes through $0$ or
% the coefficients are singular, but solutions to
% nonlinear ODE\kern .5pt s may have all kinds of singularities, which
% may be {\em movable} in the sense that they
% appear at locations dependent on initial data.
% A natural question for a mathematician is, are there
% nonlinear ODE\kern .5pt s whose movable singularities are only
% poles, never branch points?  It turns out that a
% full answer to this question is known in the case of second-order ODE\kern .5pt s
% of the form $y'' = F(t,y,y')$, where $F$ is a rational function.
% Paul Painlev\'e showed that all such equations possessing
% this property that the movable singularities are poles
% can be organized into 50 classes, 44 of which
% can be reduced to other known functions.\footnote{Painlev\'e
% was not your average mathematician.  In 1917 and again
% in 1925, he served as Prime Minister of France, and he is
% buried in the Panth\'eon in Paris.}  The remaining
% six classes are called the {\bf Painlev\'e equations,} and
% their solutions are known as {\bf Painlev\'e transcendents.}
% The {\bf Painlev\'e I} equation is 
% $$ y'' = 6y^2 + t, \eqno (21.5) $$
% and the {\bf Painlev\'e II} equation is
% $$ y'' = 2y^3 + ty + \alpha, \eqno (21.6) $$
% where $\alpha$ is a parameter.
% </latex>

%%
% <latex>
% Let us take a look at (21.5), the Painlev\'e I equation.
% In Chapter 3, we saw that
% the first-order ODE $y'=y^2$ and various generalizations
% have solutions which blow up to $\infty$.
% The blowup points are simple
% poles.  For (21.5), a second-order equation of a similar form,
% one gets double poles instead of simple poles.   For example,
% here is a solution to (21.5) on the interval $[-15,2.27\kern .5pt]$
% with ``middle conditions'' $y(0) = y'(0) = 0$.  For
% $t<0$, the solution is oscillatory and nonsingular, but
% for $t>0$ the curve is approaching a double pole at $t\approx 2.6$.
% </latex>
N = chebop(0,2.27); N.op = @(t,y) diff(y,2)-6*y^2-t;
N.lbc = [0;0]; y = N\0; plot(y,CO,ivpnl)
N = chebop(-15,0); N.op = @(t,y) diff(y,2)-6*y^2-t;
N.rbc = [0;0]; y = N\0; hold on, plot(y,CO,ivpnl), hold off
axis([-15 3 -4 10])
title('Fig.~21.8.~~Solution to Painlev\''e I eq.~(21.5)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The reader of this chapter will suspect that if $y(t)$
% oscillates like this on the negative $t$-axis, there must
% be singularities nearby in the left half of the complex
% $t$-plane.  This is true, and in fact, there is an infinite
% array of double poles extending to $\infty$ in all directions.
% </latex>

%%
% <latex>
% The Painlev\'e I equation does have smooth solutions for certain
% boundary data, which we can isolate most easily by
% solving a BVP.\ \ Here is an example:
% $$ y'' = 6y^2 + t, \quad t \in [-24,0\kern .3pt], ~ y(-24) = 2, ~y(0) = 0.
% \eqno (21.7) $$
% </latex>
N = chebop(-24,0); N.op = @(t,y) diff(y,2) - 6*y^2 -t;
N.lbc = 2; N.rbc = 0; y = N\0; plot(y,CO,bvpnl)
title('Fig.~21.9.~~Solution to Painlev\''e I eq.~(21.7)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Because of the exponential sensitivities introduced by the
% nonlinear term, one is hardly likely to find such smooth
% solutions by solving an IVP.\ \ To illustrate this
% sensitivity, here we superimpose
% on the figure the solution to 
% an IVP (or if you prefer a {\em final-value problem})
% corresponding to the same solution just plotted, except that the boundary
% conditions are both specified at the right-hand boundary, with
% a small perturbation introduced in the derivative:
% $$ y'' = 6y^2 + t, \quad t \in [-24,0\kern .3pt], ~ y(0) = 0, ~y'(0)=
% 0.999999\alpha.
% \eqno (21.8) $$
% Here $\alpha\approx -0.451427$ is the value of $y'(0)$ corresponding to the solution
% of (21.7).
% </latex>
yp = diff(y); alpha = yp(0); N.lbc = [];
N.rbc = [0; 0.999999*alpha];
y = N\0; hold on, plot(y,CO,ivpnl), hold off
title('Fig.~21.10.~~Another nearby solution superimposed',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The new solution matches the smooth one approximately for
% $t\in [-5,0\kern .3pt]$, but then diverges to an oscillatory form.
% </latex>

%%
% <latex>
% It is impossible to understand Painlev\'e equations very fully
% by looking just on the real axis, however, as is shown beautifully
% by the complex plane explorations of the paper
% by Fornberg and Weideman cited as our favorite reference below.
% </latex>

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

%%
% <latex>
% Let $m\in [\kern .3pt 0,1)$ be a parameter and consider the nonlinear
% second-order initial-value problem
% $$ y'' = -(1+m)y + 2m y^3, \quad y(0) = 0, ~ y'(0) = 1 .
% \eqno (21.9) $$
% If $m$ is zero, the equation is just $y'' = -y$ and the
% solution is $y(t) = \sin(t)$.  As~$m$ increases, the
% nonlinearity becomes stronger.  Here is what we find for
% $m = 0.998$ over the interval $t\in [\kern .3pt 0,100\kern .3pt ]$.
% </latex>
N = chebop(0,100); N.lbc = [0;1];
m = 0.998; N.op = @(y) diff(y,2) + (1+m)*y - 2*m*y^3;
y = N\0; plot(y,LW,1,CO,ivpnl), ylim([-2 2]), grid on
title(['Fig.~21.11.~~Solution of (21.9): ' ...
       'Jacobi sine function $\hbox{sn}(t,0.998)$'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Like $\sin(t)$, this function oscillates
% between $-1$ and $1$, but its shape is squarer.
% The period $T$ is about three times longer than for $\sin(t)$,
% reflecting how much time this function spends lingering
% near $\pm 1$:
% </latex>
[~,maxima] = max(y,'local'); T = maxima(3) - maxima(2)

%%
% <latex>
% Now $\sin(t)$ is completely smooth, analytic throughout
% the complex $t$-plane.  The function $y(t)$, on the other hand,
% makes rapid transitions between $-1$ and~$1$, which will become
% more abrupt as $m$ increases toward~1.\ \ We may accordingly guess that if
% this function is analytically continued into the complex plane, there will be
% singularities above and below the $t$-axis near $t=0, T/2, T, \dots .$
% A call to {\tt aaa} confirms this prediction.  In the next figure,
% the dots show the poles of the {\tt aaa} rational approximant, and
% the vertical lines mark ${\hbox{Re}\kern .3pt}(t) = 0, T, \dots , 6T$.
% </latex>
[r,pol] = aaa(y,'tol',1e-8);
plot([-10 110],[0 0],LW,.7,CO,.75*[1 1 1]), 
grid off, xlim([-10 110]), axis equal, hold on
for k = 0:6
  plot(k*T*[1 1],[-30 30],LW,.7,CO,.75*[1 1 1])
end
plot(pol,'.r',MS,10)
title(['Fig.~21.12.~~Poles near the real axis ' ...
       'found by {\tt aaa}'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% This is a striking configuration.  In fact, the function $y$ defined
% by (21.9) is {\em doubly periodic} in the complex $t$-plane,
% periodic not only in the real but also the imaginary direction.
% (The poles of the mathematically exact function
% keep going forever in a perfectly regular array,
% but {\tt aaa} has just estimated some of 
% the poles near $[\kern .3pt 0 , 100\kern .3pt]$.)
% It is the function known as the {\em Jacobi sine function} with
% parameter $m = 0.998$, written $\hbox{sn}(t,0.998)$.
% Although it is analytic for real values of $\kern .3pt t$, it has poles for
% complex~$t$, and they lie in an infinite doubly periodic array.
% Here is the period in the imaginary direction (close to $\pi$, but
% different):
% </latex>
imagpol = sort(abs(imag(pol)));
T2 = 2*min(abs(imag(pol)))

%%
% <latex>
% \noindent
% If we include horizontal lines in the plot too, we
% get an image marking fundamental domains of periodicity.
% </latex>
for k = -10:10
  plot([-10 110],k*T2*[1 1],LW,.7,CO,.75*[1 1 1])
end
plot(pol,'.r',MS,10), hold off
title(['Fig.~21.13.~~Lines mark both real and ' ...
       'imaginary periods'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Here, arbitrarily, we verify the double periodicity
% by comparing $r(t)$ for $t = 5-i$ and the same value plus
% $T$, $iT_2$, and $T+iT_2$.  (See Exercise 21.8.)
% </latex>
t = 8-1i;
r(t+[0; T; 1i*T2; T+1i*T2])

%%
% <latex>
% Doubly periodic functions, also known as {\em elliptic functions}, are
% among the core topics of complex analysis.
% </latex>

%%
% <latex>
% \smallskip
% {\sc History.}  The idea of giving physical meaning to imaginary
% space or time has proved fruitful for quite a few
% problems of engineering and physics, from the design of transonic airplane wings 
% to the explanation of the Big Bang.  Stephen Hawking presented the latter idea
% in his 1988 blockbuster book {\em A Brief History of Time}.  The elegant,
% singularity-free frame of reference in which to understand the
% universe, Hawking suggests,
% is one in which $\kern .3pt t\kern .3pt$ runs in the
% imaginary direction.  In that direction, he
% proposes, the universe's boundary conditions are periodic,
% so there are no boundaries and
% no singularities.  It's only if you turn
% a right angle and do analytic continuation into the real-$t$ direction --- 
% analytically continuing
% back to around $t = -\hbox{14,000,000,000}$ years, to be precise ---
% that you encounter the mother of all singularities.
% \smallskip
% </latex>

%%
% <latex>
% {\sc Our favorite reference.}  There is an excellent classic
% book by Einar Hille on {\em Ordinary Differential Equations in the Complex Domain,}
% but our favorite reference for this topic, with spectacular figures
% illustrating the functions in question, is Fornberg and
% Weideman, ``A numerical methodology for the Painlev\'e equations,''
% {\em Journal of Computational Physics,} 2011.
% </latex>

%%
% <latex>
% \smallskip
% \begin{displaymath}
% \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl
% {\sc Summary of Chapter 21.}  Many ODE\kern .5pt s
% make sense for complex as well as real arguments.
% Solutions\/ $y(t)$ for complex\/ $t$ can be defined by applying the
% ODE in the complex plane, or by analytic continuation from 
% solutions for real\/ $t$.  If $y(t)$ changes rapidly as\/ $t$
% varies through certain real values, there are usually one
% or more singularities of\/ $y$ nearby in the complex $t$-plane.
% Typically such singularities will be branch points, requiring
% associated branch cuts, but in the case of Painlev\'e equations,
% the only movable singularities (i.e., lying at data-dependent
% locations) are poles.
% \vspace{2pt}}}
% \end{displaymath}
% </latex>

%%
% <latex>\smallskip\small\parskip=1pt\parindent=0pt
% {\em Exercise $21.1$.  Painlev\'e I change of variables.}
% Suppose $y(t)$ is a solution to (21.5).  Show that another
% solution is $\omega^3 \kern .3pt
% y(\omega\kern .3pt t)$, where $\omega = \exp(2\pi i/5)$.
% \par 
% {\em \underline{Exercise $21.2$}.  Multiple solutions to Painlev\'e BVP.}\ \ Compute
% and plot another solution to (21.7) by starting from the
% initial guess $y(t) = -\exp(-10(t+1)^2)$.  What is its minimum value?
% \par
% {\em \underline{Exercise $21.3$}.  Some singularities are just pseudo-singularities.}
% {\em (a)} Repeat the calculation of Figure 20.4 for the ODE of (20.6),
% $\varepsilon y'' + x\kern .1pt y' + x\kern .1pt y = 0$ with $\varepsilon = 0.001$.
% Then plot the poles of the AAA approximants of the solution $y$, as in Figures
% 21.7 and 21.12, for ${\tt tol} = 10^{-4}$, $10^{-6}$, and $10^{-8}$.  
% You will see that the poles of these approximations vary with the
% tolerance.  They are not approaching any actual singularities of the
% true solution $y$, for $y$, despite
% its steep interior layer near $x=0$, cannot have any singularities since the 
% equation is linear with analytic coefficients and a nonzero
% leading-order coefficient.
% {\em (\kern .7pt b)} In such a case
% $y$, though analytic as a function of $x$, must increase
% rapidly in magnitude as $x$ leaves the real axis.  (One can prove
% this via {\em Cauchy's estimate.})  To confirm
% this, examine values of $x$ along the
% imaginary axis.  First, make the substitution $x = i s$
% as in (21.3) and write down the form that the ODE takes when
% transformed to the $s$ variable.
% {\em (c)} Now solve (20.6) numerically
% for $x \in [\kern .3pt 0, 0.2 \kern .3pt i \kern .2pt ]$, that is,
% $s \in [\kern .3pt 0, 0.2 \kern .2pt ]$,
% taking as initial data the appropriately
% transformed values $y(0)$ and $y'(0)$ from
% the solution of {\em (a)}.  Make a semilogy plot of
% $|y(x)|$ against $|x|$.  How big is $|y(0.2\kern .3pt i)|\kern .5pt ?$
% \par
% {\em \underline{Exercise $21.4$}.  Bypassing blowup.}  As we have discussed
% at several points in the book, positive 
% solutions to $y' = y^2$ blow up to $\infty$ in finite time.  The situation
% changes, however, if $y$ is even very slightly complex.
% {\em (a)} Draw a quiver plot of $y' = y^2$ in the
% complex $y$-plane and describe the shape of the complex
% trajectories in this plot.
% {\em (\kern .7pt b)} Solve $y' = y^2 + f$ with $y(0)=1$
% for $t\in [\kern .3pt 0,1]$, where
% $f$ is the complex random function generated by {\tt rng(1)},
% \verb| f = 0.01*randnfun([0 1],'norm','complex')|.
% Plot the solution with {\tt axis equal}.  (For
% analysis of such problems see
% Herzog and Mattingly, ``Noise-induced stabilization of
% planar flows I'', {\em Electronic Journal of Probability,} 2015.)
% \par
% {\em \underline{Exercise $21.5$}.  Jacobi sine function.}
% A chebfun {\tt y2} representing the Jacobi sine function
% as plotted in Fig.~21.11 could be constructed directly
% from the Matlab function \verb|ellipj(t,0.998)|.
% Do this and plot {\tt abs(y-y2)} to give an indication of accuracy of
% this Chebfun ODE solution.
% </latex>
