%% Appendix A. Chebfun and its ODE algorithms

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

%%
% This is not a book of numerical analysis, yet its
% heart is numerical explorations.
% Here we outline Chebfun's ODE algorithms, with
% references to papers giving further details listed at
% the end, and mention some of the people who developed them.

%%
% <latex>
% \smallskip
% {\em Functions.}
% The Chebfun project began at Oxford in 2002 with Trefethen and his 
% DPhil student Zachary Battles [3,20].  The starting idea, which
% remains the project's central vision, is to utilize 
% Chebyshev interpolants and expansions to enable
% ``numerical computing with functions.''  Specifically, Chebfun
% overloads Matlab commands for discrete vectors 
% to analogues for continuous functions, aiming to deliver speedy results 
% that are accurate to close to machine precision.
% </latex>

%%
% In Chebfun, each function (or piece of a function, for
% piecewise representations) is represented by a polynomial whose
% domain is a real interval, $[-1,1]$ by default.  Polynomials
% are manipulated via
% Chebyshev series and coefficients defined by the formulas
% $$ f(x) = \sum_{k=0}^\infty a_k T_k(x), \quad
% a_k = {2\over \pi} \int_{-1}^1 {f(x) T_k(x)\over \sqrt{1-x^2}}\kern 
% 1pt dx, \eqno ({\rm A}.1) $$
% where $T_k$ is the degree $k$ Chebyshev polynomial (for
% $k=0,$ the factor $2/\pi$ changes to $1/\pi$).
% If $f$ is Lipschitz continuous, the series converges uniformly and
% absolutely, and the smoother $f$ is, the faster the convergence.
%%
% For example, $f(x) = 50x\exp(-200x^2)-\tanh(e^x-1)$ 
% can be constructed and plotted like this:
ODEformats
f = chebfun(@(x) 50*x*exp(-200*x^2)-tanh(exp(x)-1));
plot(f), ylim([-1.8 1.8]), n = length(f)-1;
title(['Fig.~A.1.~~A chebfun $f$ --- in this ' ...
       'case a polynomial of degree ' int2str(n)],FS,11)
%%
% <latex> \vskip 1.02em </latex>
%% 
% <latex>
% \noindent
% The approximation {\tt f}
% is called a chebfun, with a lower-case c.
% To compute quantities like the maximum, the integral, or the norm,
% Chebfun uses further Chebyshev-based algorithms,
% all invoked by overloads of familiar Matlab commands.
% </latex>
format long, max(f), sum(f), norm(f)
%%
% <latex>
% \noindent
% Algebraic and other operations are also
% overloaded in the expected manner: 
% </latex>
min(exp(sin(1/(1+f^2))))

%%
% <latex>
% \noindent
% Zeros of functions are computed by the command {\tt roots:}
% </latex>
roots(f)

%%
% <latex>
% \noindent
% The relative accuracy of each computation is usually about
% 16 digits, and in principle,
% the user need have no knowledge of the underlying algorithms.  For
% example, this function
% $f$ is approximated by a polynomial of degree in the hundreds.
% Here are
% the absolute values of its Chebyshev series coefficients, plotted
% on log scale.  (As it happens, for this function the even coefficients
% decrease much faster than the odd ones, eventually reaching a 
% plateau around $10^{-16}$ because of rounding errors.)
% </latex>
plotcoeffs(f,CO,'k',MS,6), ylim([1e-20 1e2])
title(['Fig.~A.2.~~Absolute values of ' ...
       'Chebyshev coefficients of $f$'],FS,11)
xlabel('$k$',FS,10), ylabel('$|a_k|~~~~$',FS,10,'rotation',0)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% To find this approximation, Chebfun samples $f$ on
% Chebyshev grids $\{x_j^{(n)}\}$ with
% $n+1 = 17, 33, 65, \dots$ points, defined by
% $$ x_j^{(n)} = -\cos(\kern .5pt j\pi/n),\quad 0\le j \le n. \eqno ({\rm A}.2) $$
% On each grid, it
% determines the Chebyshev series coefficients of the degree $n$ polynomial
% interpolant through the samples.  When the coefficients
% hit a plateau of rounding errors at a relative magnitude
% of about $10^{-16}$, the grid refinement stops and
% the series is trimmed~[1].  For this function, the
% plateau is first detected on the grid of 257 points, and we can
% see what the Chebyshev series looks like before trimming
% by instructing Chebfun to sample in exactly 257 points: 
% </latex>
f = chebfun(@(x) 50*x*exp(-200*x^2)-tanh(exp(x)-1),257);
plotcoeffs(f,CO,'k',MS,6), axis([0 270 1e-20 1e2])
title('Fig.~A.3.~~ Chebyshev coefficients on a 257-point grid',FS,11)
xlabel('$k$',FS,10), ylabel('$|a_k|~~~~$',FS,10,'rotation',0)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% The mathematics behind most Chebfun algorithms is presented
% in the book {\em Approximation Theory and Approximation
% Practice} [19], and user information can be found in
% the Chebfun {\em Guide} [10] and web site [6].
% The collection of hundreds of Examples posted at [6] may be
% particularly useful.
% </latex>

%%
% <latex>
% Two new team members joined the Chebfun project
% during 2006--07, Ricardo Pach\'on and Rodrigo Platte, and they
% introduced piecewise representations for problems
% with discontinuities [14], which were later
% incorporated in the BVP solution algorithms by Nick Hale
% and Toby Driscoll.
% Piecewise representations are exploited at many points 
% in this book, starting with the problems
% (2.15) (``Sydney Opera House'') and (2.16) (``Batman'').
% </latex>

%%
% <latex>
% \smallskip
% {\em Linear BVPs.}
% The same Chebyshev grids and expansions that work so well for
% approximating functions also have a distinguished history
% for solving differential equations, where they go
% by the name of {\em Chebyshev spectral methods} [18].
% Beginning in 2008, following an initial proposal by Folkmar
% Bornemann, Toby Driscoll created the BVP side of Chebfun, initially
% for linear problems [8].  The principle is the same as before:
% solve the problem on grids of size 17, 33, 65 (actually, slightly
% different parameters are used), and in each case, examine the Chebyshev
% series for convergence.
% When a grid is found for which convergence is achieved, the
% series is trimmed as usual, and the corresponding
% chebfun is returned as a solution to the BVP.\ \ Each solution
% on a given grid involves
% a discretization by so-called Chebyshev discretization
% matrices.  Initially, Chebfun followed the traditional
% ``square matrix'' discretization strategy described
% in [18] (and many other places),
% but later, Driscoll and Hale found that a
% different ``rectangular matrix'' formulation offers greater
% reliability and robustness, especially for problems with
% nonstandard boundary conditions or multiple dependent
% variables [2,9,22].
% </latex>

%%
% <latex>
% For example, here is the solution to a ``wavy Airy function''
% problem adapted from Chapter~7:
% </latex>
L = chebop(-50,20);                            % domain
L.op = @(x,y) diff(y,2) - (2+sin(x/2))*x*y;    % diff'l operator
L.lbc = 0; L.rbc = 0;                          % boundary conds
y = L\1;                                       % solution
plot(y,LW,.6,CO,bvp)                           % plot
title('Fig.~A.4.~~A wavy Airy function',FS,11)
%%
% <latex> \vskip 1.02em </latex>
%%
% <latex>
% \noindent
% This polynomial representation is of degree over 500:
% </latex>
plotcoeffs(y,CO,'k',MS,5)
title('Fig.~A.5.~~Chebyshev coefficients of a wavy Airy function',FS,11)
xlabel('$k$',FS,10), ylabel('$|a_k|~~~~$',FS,10,'rotation',0)
%%
% <latex> \vskip 1.02em </latex>
%%
% <latex>
% \noindent
% Since the problem domain is $[-50,20\kern .3pt]$ rather than $[-1,1]$,
% the coefficients plotted correspond to
% Chebyshev polynomials appropriately transplanted.
% </latex>

%%
% <latex>
% In the code segment above, the command that invokes
% the solution of the BVP is \verb|y = L\1|: Matlab backslash.
% This is a continuous analogue of Matlab's use of backslash
% to solve a matrix problem ${\bf A}{\bf x}={\bf b}$ via \verb|x = A\b|.  In both
% cases, the purpose of the compact notation is to highlight the
% conceptually simple solve operation required while suppressing
% algorithmic details.
% Just as an advanced Matlab user
% can call {\tt linsolve} instead of \verb|\| to specify nondefault
% algorithmic parameters, an advanced Chebfun user can call
% {\tt solvebvp} instead of~\kern .5pt\verb|\|\kern .5pt, though we have not found
% the need for that in this book.
% </latex>

%%
% An interesting conceptual step was involved in advancing
% Chebfun from functions to solutions of BVPs.  To represent a
% function, one simply samples it at more and more points.
% For BVPs, however, the function being sampled is not known
% a priori.
% The ``samples'' at each particular point
% will vary from grid to grid until the spectral
% approximation converges, and moreover, the work involved
% on each grid
% is proportional to the cube of the number of grid points.
% These new features of the problem have certain engineering
% consequences but in the end are not so significant.

%%
% <latex>
% \smallskip
% {\em Eigenvalue problems.}
% Besides linear BVPs, in 2008 Driscoll also introduced Chebfun commands
% {\tt eigs} for eigenvalue problems (Chapter 6), {\tt expm}
% for exponentials of linear operators,
% and {\tt fred} and {\tt volt} for Fredholm and Volterra
% integral equations [7].  Like \kern .5pt\verb|\|\kern .5pt,
% these operations are implemented via Chebyshev discretization matrices of
% adaptively determined dimensions.
% Here, to illustrate (compare Fig.~6.8),
% are the first six eigenfunctions of
% the harmonic oscillator of quantum mechanics, the Schr\"odinger
% equation $-u'' + x^2 u = \lambda u$ on $[-5,5]$ with 
% $u(\pm 5) = 0$ (an approximation to the infinite real line):
% </latex>
L = chebop(-5,5); L.op = @(x,u) -diff(u,2) + x^2*u;
L.lbc = 0; L.rbc = 0; [V,D] = eigs(L,6); plot(V)
title('Fig.~A.6.~~Eigenfunctions of the harmonic oscillator',FS,11)
%%
% <latex> \vskip 1.02em </latex>
%%
%%
% <latex>
% \noindent
% In standard Matlab, users of {\tt eigs} can specify
% whether they want to compute eigenvalues of largest magnitude
% (the default),
% largest real part, largest imaginary part, etc.\ \ In Chebfun
% the same options are available, but the default is a
% different choice entirely: eigenvalues associated
% with {\em smoothest eigenfunctions.}  In applications, like
% the one just illustrated, these are typically the
% eigenfunctions of interest, and this small point of eigenvalue
% calculations illustrates how new considerations
% may come into play when one takes the step from discrete to
% continuous.
% </latex>

%%
% <latex>
% \smallskip
% {\em Nonlinear BVPs.}
% In 2009 two further Chebfun graduate students
% joined the project, Nick Hale and
% \'Asgeir Birkisson.
% Birkisson began to work
% with Driscoll on the challenge of extending Chebfun's linear BVP
% algorithms to nonlinear problems.  Nonlinearity requires a Newton iteration,
% as mentioned in Exercise~5.10 and Chapter 16,
% and in the Chebfun spirit, we wanted to realize this in
% ``continuous mode'' --- or as it is often put in computational
% science, in the mode of {\em solve-then-discretize\/}
% rather than {\em discretize-then-solve.} 
% The reader will be aware that Newton's method for finding
% a root of a nonlinear equation $f(y)$ requires evaluation of the
% derivative $f'(y)$ at various points.  Similarly, for a
% system of equations ${\bf f}({\bf y}) = {\bf 0}$,
% one needs to evaluate a
% Jacobian matrix of partial derivatives, $J_{ij} = \partial f_i/\partial 
% y_j$.
% In Chebfun's continuous setting,
% this matrix becomes an infinite-dimensional Fr\'echet
% derivative linear operator
% that must be discretized.  Birkisson and Driscoll found a
% way to achieve this using automatic
% differentiation [5], and together with Nick Hale, they became
% the primary developers of the ODE side of Chebfun for the
% next five years.  The linear part of BVP solution is
% now contained in the Chebfun class {\tt linop,} which the
% user normally does not call directly.
% </latex>

%%
% <latex>
% Even for nonlinear problems, Chebfun's basic BVP command is
% backslash.  One can use {\tt solvebvp,} however, to output
% additional information about the solution.
% For example, here we use {\tt solvebvp} to solve
% the problem (16.8),
% </latex>
N = chebop(-1,1); N.lbc = 0; N.rbc = 0; 
N.op = @(x,y) 0.2*diff(y,2) + y + y^2; 
N.init = chebfun('sin(pi*x)');
[y,info] = solvebvp(N,1); plot(y,CO,bvpnl), ylim([-1.6 1.6])
title('Fig.~A.7.~~Solution to (16.8)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% One of the fields of {\tt info} then enables us to track
% the norms of the Newton updates during the iteration, showing
% 8 steps and quadratic convergence.
% </latex>
semilogy(info.normDelta,'.-k'), ylim([1e-18 1e2])
title('Fig.~A.8.~~Norms of Newton updates during iteration',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \smallskip
% {\em Initial-value problems.}
% By 2013, Chebfun was a powerful and user-friendly tool
% for solving ODE BVPs and associated eigenvalue problems,
% and besides simplicity, three of its
% impressive features, though we have not emphasized
% them in this discussion, were its treatment of systems
% of equations, problems with discontinuities, and nonstandard
% ``boundary'' conditions such as integral conditions.  Something
% big was still missing before this book could be written, however,
% and that was an equally effective treatment of IVPs following
% the same syntax.
% </latex>

%%
% <latex>
% Algorithmically, IVPs and BVPs differ greatly, for the best
% general methods for solving IVPs, whose roots
% go back to Adams, Bashforth, Runge, and
% Kutta more than a century ago, are based on {\em marching}
% in time, in contrast to the global spectral methods
% associated with Chebfun backslash.  
% Global methods can be used for IVPs, and for linear problems
% this is typically quite effective, but for nonlinear problems
% the introduction of a Newton iteration is terribly wasteful
% and generally leads to nonconvergence or at best to great slowdown.
% Using backslash, Chebfun would have no chance of producing the
% solutions of the Lorenz and van der Pol equations shown at many
% points in this book, for example.
% (Conversely, using a shooting iteration based on
% marching as mentioned in Chapter~16,
% one would have trouble with many BVPs.)
% </latex>

%%
% <latex>
% In Matlab, the disjunction between BVPs and IVPs is reflected
% in the provision of the separate codes {\tt bvp4c, bvp5c}
% for the former and
% {\tt ode45, ode113, ode15s}, etc.\ for the latter.   These
% use different algorithms and different syntaxes.  For Chebfun,
% Platte and Hale introduced Chebfun overloads
% of {\tt ode45} and {\tt ode113}, among others, in 2008--09.  The most important
% Matlab ODE IVP solver for
% our purposes is {\tt ode113,} because it is usually most
% efficient for the high-accuracy solutions aimed for in Chebfun (except
% for very stiff problems, where {\tt ode15s} is superior; for an example
% see the Oregonator in Appendix~B).
% This code is based on an Adams linear multistep
% method with variable order ranging from $1$ to $13$.
% Hale's algorithmic idea was to solve an IVP by calling Matlab's
% {\tt ode113} with a tight tolerance specification, then convert
% the result to a chebfun by the usual process of sampling on finer
% and finer grids.
% </latex>

%%
% <latex>
% For many purposes, however, including a book like this, it would not
% do to confuse users with
% one syntax for BVPs and another for IVPs.
% We had to unify the
% disparate algorithms within the backslash framework,
% and here, a seemingly small obstacle proved substantial:
% whereas Chebfun and its spectral methods readily discretize
% BVPs of any order, {\tt ode113} expects every IVP to be
% formulated as a first-order system.  So a method was needed
% to convert higher-order ODE specifications automatically
% to first-order form, and in 2014, Birkisson achieved this
% by a method based on operator overloading and the generation
% of syntax trees of expressions as mathematical programs are
% evaluated [4].
% </latex>

%%
% For example, here is the solution of
% a van der Pol equation with a time-varying coefficient:
N = chebop(0,40);
N.op = @(t,y) diff(y,2) - 0.2*t*(1-y^2)*diff(y) + y;
N.lbc = [3; 0]; y = N\0; plot(y,CO,ivpnl)
title('Fig.~A.9.~~Time-dependent van der Pol IVP',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The solution is a single long chebfun, a polynomial of
% high degree:
% </latex>
plotcoeffs(y,CO,'k',MS,3)
xlabel('$k$',FS,10), ylabel('$|a_k|~~~~$',FS,10,'rotation',0)
title(['Fig.~A.10.~~Chebyshev coeffs for ' ...
       'time-dependent van der Pol IVP'],FS,11)
format short
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Thus finally, beginning in 2014, a Chebfun
% user could specify a BVP or an IVP using
% the same backslash syntax and get a
% high-accuracy result quickly. The algorithms are completely
% different in the two cases, but the inputs have
% the same syntax and the outputs are always chebfuns.
% Advanced users can fine-tune the computation with
% {\tt solvebvp} or {\tt solveivp}, and an IVP can be solved
% by global spectral methods rather than marching, if desired,
% by calling {\tt solvebvp}.
% </latex>


%%
% <latex>
% \smallskip
% {\em Periodic problems.}
% Another new feature also appeared in 2014, introduced by
% new Chebfun contributor Grady Wright, who was visiting Oxford on sabbatical
% from Boise State University.
% From the beginning, Chebfun had been based on nonperiodic Chebyshev
% discretizations, even though one could specify a \verb|'periodic'|
% flag as a boundary condition for a BVP.\ \ This
% was not really a satisfactory way
% to deal with periodic problems, however, since Chebyshev matrices have
% large condition numbers associated with the clustered
% grids, and derivatives of Chebyshev representations
% of periodic functions quickly develop discontinuities across
% the boundary.  In 2014, Wright introduced Fourier-based ``trigfun''
% representations to Chebfun and the corresponding Fourier
% as opposed to Chebyshev discretizations for BVPs [21].  Together
% with Oxford graduate student
% Hadrien Montanelli, Wright adapted Chebfun's Chebyshev-based
% BVP capabilities to Fourier analogues for periodic functions,
% as presented in Chapter~19.
% </latex>

%%
% <latex>
% \smallskip
% {\em PDEs in space and time.}
% As illustrated in Chapter 22, besides ODE\kern .5pt s, Chebfun can also solve PDEs
% involving time $t$ and one space variable $x$.
% The general tool for such solutions, {\tt pde15s}, was developed by Nick
% Hale beginning in 2009, with further improvements later
% by \'Asgeir Birkisson.  This code uses spectral discretization in
% $x$ combined with Matlab's time-stepper {\tt ode15s} in $t$, 
% adapting grids and refining time steps adaptively to achieve 
% chebfun outputs.  Both nonperiodic
% (Chebyshev) and periodic (Fourier) discretizations are
% available.  Unfortunately, there
% is no publication describing {\tt pde15s}.
% </latex>

%%
% <latex>
% For time-dependent PDEs on a periodic 1D space domain,
% another more specialized option
% was added by Hadrien Montanelli in 2015 [12]:
% the code {\tt spin,} which stands
% for {\em stiff PDE integrator.}  This code uses exponential integrator
% formulas to solve equations of the form $u_t^{} = Lu + N(u)$, where
% $L$ is a linear differential operator and $N$ is a nonlinear differential
% or algebraic operator of lower order.  Equations of this type include
% the Burgers, Korteweg--de Vries, nonlinear Schr\"odinger,
% FitzHugh--Nagumo, Allen--Cahn, Cahn--Hilliard, Gray--Scott,
% Nikolaevskiy, and Kuramoto--Sivashinsky equations, and
% demos for the examples just listed are available with a syntax
% like \verb|spin('kdv')|.
% </latex>

%%
% <latex>
% \smallskip
% {\em Chebgui.}
% An easy way for users to explore all of these many sides of
% Chebfun's differential equations capabilities is with the
% graphical user interface {\tt chebgui,} written by
% Birkisson and initially also Hale.  In Chebfun, just type
% {\tt chebgui} to get started, and note the many example problems
% available under the Demos tab.  Chebfun code can be generated
% with the ``Export to m-file'' button, a good starting point
% for more finely-tuned computations.
% </latex>

%%
% <latex>
% \smallskip
% {\em Additional capabilities and higher dimensions.}
% This completes our outline of the 
% solution of ODE\kern .5pt s in Chebfun.  We have given little attention to
% an aspect of the problem that in practice adds a great deal of
% complexity to the design: treatment of systems of equations, with
% its associated chebmatrix class in the software.
% We have also not mentioned various additional features
% of Chebfun
% that are numerically very interesting, though they have
% not played a role in this book, including first-kind as opposed
% to second-kind Chebyshev grids, ultraspherical discretizations,
% ODE\kern .5pt s with unknown parameters, differential-algebraic equations
% (DAEs), and unbounded domains.
% Chebfun has capabilities in all of these areas.
% </latex>

%%
% <latex>
% There is also a big part of Chebfun that does not relate directly
% to this book: two- and three-dimensional spatial domains.
% Such capabilities began with Chebfun2, designed by graduate
% student Alex Townsend around 2011--12 [15], and a 3D extension
% Chebfun3 was released by postdoc Behnam Hashemi in 2016 [11].
% Building on the periodic trigfun representations, a new
% Spherefun class for computing on spheres was introduced in 2015
% by Townsend, Wright, and Heather Wilber, a student of Wright's
% at Boise State [16].  Later at Cornell, 
% Wilber went on to produce a related
% Diskfun class in 2016 for computing on a disk, again in collaboration
% with Townsend and Wright [17].  All four of these
% multidimensional parts of Chebfun --- Chebfun2, Chebfun3, 
% Spherefun, and Diskfun --- come with PDE capabilities, including
% solution of reaction-diffusion and other stiff PDEs with
% {\tt spin2}, {\tt spin3}, and {\tt spinsphere} [13].
% </latex>

%%
% <latex>
% \smallskip
% {\em Numerical limitations.}
% Finally, we must emphasize that Chebfun has limitations.
% When the ODE side of Chebfun was developed, beginning in 2008,
% the expectation was not that it would compete with
% existing software, merely that it would
% be a convenient tool for exploring certain algorithms
% in a continuous setting.  We were surprised
% to find how useful the tool was in practice and to find
% that in some cases it was competitive after all,
% especially for BVPs solved to high accuracy.
% Still, by pushing Chebfun too far, it is easy to make it slow down,
% break, or at least require expert handling, as we now
% outline following the paragraph at the end of Chapter~1.
% {\em Stiffness.}  If an ODE features a wide range of time
% scales, Chebfun is likely to have trouble.
% For example, the computation will become very slow if
% the coefficient of the van der Pol equation (1.2) is changed from
% $0.3$ to $0.01$.  For such problems one can switch from the
% default {\tt ode113} to the stiff solver {\tt ode15s}, as illustrated
% in equation~36 in Appendix~B.\ \ {\em Scaling.}  Some
% ODE\kern .5pt s feature variables that differ by 
% many orders of magnitude from $1$, from other variables, or
% from values taken by the same variable at other times. 
% Taking advantage of floating
% point arithmetic, some numerical ODE software can cope with such problems,
% but Chebfun will have difficulty in its default mode.
% For a simple problem
% like $y' = -y$ on $[\kern .3pt 0,1]$ with initial condition
% $y(0) = 10^{-20}$, one can get a successful result by overriding
% the default absolute tolerance, as illustrated in
% equation~xx of Appendix~B.\ \ For a
% more difficult case like the same equation on $[\kern .3pt 0,
% 50]$, Chebfun will be sure to fail because its function representations
% are global.
% {\em Discontinuities in $y$.}  Chebfun can treat coefficients with
% discontinuities in the independent variable, as a number of our 
% examples beginning with Fig.~2.3 have shown, but it cannot handle
% discontinuities in the dependent variable.  
% Calculating the ``bounce pass'' trajectories
% of Exercise 16.2, for example, is beyond Chebfun.
% {\em Singularities.}  Many ODE\kern .5pt s have singular coefficients or leading-order
% coefficients that pass through zero, starting with the Bessel equation
% with its zero coefficients at $x=0$, and this is a familiar chapter
% of the classical theory of ODE\kern .5pt s.  Sometimes Chebfun can handle such
% problems, but not always, and if the solution itself has singularities,
% there will almost certainly be difficulty.
% {\em Side conditions.} Some ODE software can impose side conditions such
% as a nonnegativity constraint on a dependent variable, but Chebfun cannot.
% For example, it has no very satisfactory way of
% treating the ``leaky bucket'' problem $y' = y^{1/2}$
% (see Equation~xx of Appendix~B).\ \ This limitation is
% related to the previous two.
% {\em Large systems.}  Chebfun can solve a system of
% ODE\kern .5pt s involving explicitly named
% variables such as $u,v,w$, and there is also an option involving
% ``chebmatrices'' to use indexed variables such as
% \verb|u{1},u{2},u{3}|, illustrated in Equation xx of
% Appendix~B.  However, its capabilities for working with systems
% are not very developed.  This has constrained our choice
% of illustrations
% in this book, where no systems of more than three variables are
% to be found, and it would make Chebfun an unattractive choice
% for many larger scale problems of computational science.
% </latex>

%%
% <latex>
% It would hardly do to end on a negative note, however.
% The fact is that Chebfun's efficiency is often very good, and
% despite the limitations just outlined, it is
% probably unrivalled in {\em convenience.}
% For more examples, turn the page to Appendix~B.
% </latex>

%%
% <latex>
% \bigskip
% \parindent=0pt
% {\bf References}
% \par\medskip
% {\small
% [1] J. L. Aurentz and L. N. Trefethen, Chopping a Chebyshev series,
% {\em ACM Trans.\ Math.\ Softw.} 43 (2017), 33.
% \par\vskip 2pt
% [2] J. L. Aurentz and L. N. Trefethen, Block operators and
% spectral discretizations, {\em SIAM Rev.} 59 (2017), pp.~423--446.
% \par\vskip 2pt
% [3] Z. Battles and L. N. Trefethen, An extension of MATLAB to
% continuous functions and operators, {\em SIAM J. Sci.\ Comp.}
% 25 (2004), pp.~1743--1770.
% \par\vskip 2pt
% [4] A. Birkisson, Automatic reformulation of ODE\kern .5pt s to
% systems of first order equations, to appear.
% \par\vskip 2pt
% [5] A. Birkisson and T. A. Driscoll, Automatic Fr\'echet
% differentiation for the numerical solution of
% boundary-value problems, {\em ACM Trans.\ Math.\ Softw.}\ 38
% (2012), 26.
% \par\vskip 2pt
% [6] Chebfun web site, {\tt www.chebfun.org.}
% \par\vskip 2pt
% [7] T. A. Driscoll, Automatic spectral collocation for integral,
% integro-differential, and integrally reformulated
% differential equations, {\em J. Comp.\ Phys.} 229 (2010),
% pp.~5980--5998.
% \par\vskip 2pt
% [8] T. A. Driscoll, F. Bornemann, and L. N. Trefethen,
% The chebop system for automatic solution of differential equations,
% {\em BIT Numer.\ Math.} 48 (2008), pp.~701--723.
% \par\vskip 2pt
% [9] T. A. Driscoll and N. Hale, Rectangular spectral
% collocation, {\em IMA J. Numer.\ Anal.} 36 (2016), pp.~108--132.
% \par\vskip 2pt
% [10] T. A. Driscoll, N. Hale, and L. N. Trefethen, eds.,
% {\em Chebfun Guide,} Pafnuty Publications, Oxford, 2014;
% freely available at [6].
% \par\vskip 2pt
% [11] B. Hashemi and L. N. Trefethen,
% Chebfun in three dimensions, {\em SIAM J. Sci.\ Comp.,} to appear.
% \par\vskip 2pt
% [12] H. Montanelli and N. Bootland, Solving periodic 
% semilinear stiff PDEs in
% 1D, 2D and 3D with exponential integrators, to appear.
% \par\vskip 2pt
% [13] H. Montanelli and Y. Nakatsukasa, Fourth-order
% time-stepping for stiff PDEs on the sphere,
% {\em SIAM J. Sci.\ Comp.,} to appear.
% \par\vskip 2pt
% [14] R. Pach\'on, R. B. Platte, and L. N. Trefethen,
% Piecewise-smooth chebfuns, {\em IMA J. Numer.\ Anal.} 30 (2010),
% pp.~898--916.
% \par\vskip 2pt
% [15] A. Townsend and L. N. Trefethen,
% An extension of Chebfun to two dimensions, {\em SIAM J.
% Sci.\ Comp.} 35 (2013), pp.~C495--C518.
% \par\vskip 2pt
% [16] A. Townsend, H. Wilber, and G. B. Wright, 
% Computing with functions in spherical and polar geometries I. The sphere,
% {\em SIAM J. Sci.\ Comp.} 38 (2016), pp.~C403--C425.
% \par\vskip 2pt
% [17] A. Townsend, H. Wilber, and G. B. Wright, 
% Computing with functions in spherical and polar geometries II. The disk,
% {\em SIAM J. Sci.\ Comp.,} to appear.
% \par\vskip 2pt
% [18] L. N. Trefethen, {\em Spectral Methods in MATLAB,}
% SIAM, 2000.
% \par\vskip 2pt
% [19] L. N. Trefethen, {\em Approximation Theory and Approximation
% Practice,} SIAM, 2013.
% \par\vskip 2pt
% [20] L. N. Trefethen, Computing numerically with functions instead
% of numbers, {\em Math.\ Comput.\ Sci.} 1 (2007), 9--19; expanded
% version in {\em Commun.\ ACM} 10 (2015), pp.~91--97.
% \par\vskip 2pt
% [21] G. B. Wright, M. Javed, H. Montanelli, and L. N. Trefethen,
% Extension of Chebfun to periodic functions, {\em SIAM J.
% Sci.\ Comp.} 37 (2015), pp.~C554--C573.
% \par\vskip 2pt
% [22] K. Xu and N. Hale, Explicit construction of
% rectangular differentiation matrices,
% {\em IMA J. Numer.\ Anal.} 36 (2016), pp.~618--632.
% \par}
% </latex>

