%% Appendix A. Chebfun and its ODE algorithms %% % % \setcounter{page}{293} % %% % 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. %% % % \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. % %% % 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) %% % \vskip 1.02em %% % % \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. % format long, max(f), sum(f), norm(f) %% % % \noindent % Algebraic and other operations are also % overloaded in the expected manner: % min(exp(sin(1/(1+f^2)))) %% % % \noindent % Zeros of functions are computed by the command {\tt roots:} % roots(f) %% % % \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.) % 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) %% % \vskip 1.02em %% % % \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~. 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: % 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) %% % \vskip 1.02em %% % % The mathematics behind most Chebfun algorithms is presented % in the book {\em Approximation Theory and Approximation % Practice} , and user information can be found in % the Chebfun {\em Guide}  and web site . % The collection of hundreds of Examples posted at  may be % particularly useful. % %% % % 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 , 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''). % %% % % \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} . % Beginning in 2008, following an initial proposal by Folkmar % Bornemann, Toby Driscoll created the BVP side of Chebfun, initially % for linear problems . 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  (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]. % %% % % For example, here is the solution to a wavy Airy function'' % problem adapted from Chapter~7: % 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) %% % \vskip 1.02em %% % % \noindent % This polynomial representation is of degree over 500: % 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) %% % \vskip 1.02em %% % % \noindent % Since the problem domain is $[-50,20\kern .3pt]$ rather than $[-1,1]$, % the coefficients plotted correspond to % Chebyshev polynomials appropriately transplanted. % %% % % 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. % %% % 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. %% % % \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 . 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): % 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) %% % \vskip 1.02em %% %% % % \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. % %% % % \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 , 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. % %% % % 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), % 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) %% % \vskip 1.02em %% % % \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. % semilogy(info.normDelta,'.-k'), ylim([1e-18 1e2]) title('Fig.~A.8.~~Norms of Newton updates during iteration',FS,11) %% % \vskip 1.02em %% % % \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. % %% % % 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.) % %% % % 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. % %% % % 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 . % %% % 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) %% % \vskip 1.02em %% % % \noindent % The solution is a single long chebfun, a polynomial of % high degree: % 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 %% % \vskip 1.02em %% % % \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}. % %% % % \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 . 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. % %% % % \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}. % %% % % For time-dependent PDEs on a periodic 1D space domain, % another more specialized option % was added by Hadrien Montanelli in 2015 : % 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')|. % %% % % \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. % %% % % \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. % %% % % 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 , and a 3D extension % Chebfun3 was released by postdoc Behnam Hashemi in 2016 . % 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 . 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 . 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} . % %% % % \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. % %% % % 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. % %% % % \bigskip % \parindent=0pt % {\bf References} % \par\medskip % {\small %  J. L. Aurentz and L. N. Trefethen, Chopping a Chebyshev series, % {\em ACM Trans.\ Math.\ Softw.} 43 (2017), 33. % \par\vskip 2pt %  J. L. Aurentz and L. N. Trefethen, Block operators and % spectral discretizations, {\em SIAM Rev.} 59 (2017), pp.~423--446. % \par\vskip 2pt %  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 %  A. Birkisson, Automatic reformulation of ODE\kern .5pt s to % systems of first order equations, to appear. % \par\vskip 2pt %  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 %  Chebfun web site, {\tt www.chebfun.org.} % \par\vskip 2pt %  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 %  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 %  T. A. Driscoll and N. Hale, Rectangular spectral % collocation, {\em IMA J. Numer.\ Anal.} 36 (2016), pp.~108--132. % \par\vskip 2pt %  T. A. Driscoll, N. Hale, and L. N. Trefethen, eds., % {\em Chebfun Guide,} Pafnuty Publications, Oxford, 2014; % freely available at . % \par\vskip 2pt %  B. Hashemi and L. N. Trefethen, % Chebfun in three dimensions, {\em SIAM J. Sci.\ Comp.,} to appear. % \par\vskip 2pt %  H. Montanelli and N. Bootland, Solving periodic % semilinear stiff PDEs in % 1D, 2D and 3D with exponential integrators, to appear. % \par\vskip 2pt %  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 %  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 %  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 %  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 %  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 %  L. N. Trefethen, {\em Spectral Methods in MATLAB,} % SIAM, 2000. % \par\vskip 2pt %  L. N. Trefethen, {\em Approximation Theory and Approximation % Practice,} SIAM, 2013. % \par\vskip 2pt %  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 %  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 %  K. Xu and N. Hale, Explicit construction of % rectangular differentiation matrices, % {\em IMA J. Numer.\ Anal.} 36 (2016), pp.~618--632. % \par} %