Chapter 7 described Chebfun's chebop capabilities for solving linear ordinary differential equations by the backslash command. We will now describe extensions of chebops to nonlinear problems, as well as special methods used for ODE initial-value problems (IVPs) as opposed to boundary-value problems (BVPs). Most of the design and implementation of these features was done by Asgeir Birkisson in collaboration with Toby Driscoll. For many examples and a variety of further information, see the Chebfun-based ODE textbook [Trefethen, Birkisson & Driscoll 2018], which is freely available at

10.1 Boundary-value problems: \ and solvebvp

Chebfun contains overloads bvp4c and bvp5c of MATLAB codes of the same names. However, these are not our recommended methods for solving BVPs, and we will not discuss them here. Instead, we present methods based on \ and its equivalent command solvebvp.

Recall that in Chapter 7, we realized linear operators as chebops constructed by commands like these:

L = chebop(-1, 1);
L.op = @(x,u) 0.0001*diff(u,2) + x*u;

Using such an object we can solve the BVP $$ 0.0001 u'' + xu = e^x, \qquad u(-1) = 0, ~~ u(1) = 1 $$ as follows:

L.lbc = 0; L.rbc = 1;
x = chebfun('x');
u = L\exp(x);
plot(u), grid on, ylim([-50 50])

What's going on in such a calculation is that L is a prescription for constructing matrices of arbitrary dimensions which are Chebyshev spectral approximations to the differential operator. When backslash is executed, the problem is solved on successively finer grids until convergence is achieved.

The object L is a chebop:

L =
   Linear operator:
      u |--> 0.0001.*diff(u,2)+x.*u
   operating on chebfun objects defined on:
    left boundary condition(s):
      u = 0
    right boundary condition(s):
      u = 1

Notice that Chebfun has detected that the chebop is linear. Doing this automatically is not a triviality! --- see [Birkisson & Driscoll 2013].

The same approach also works for nonlinear problems. For example, in Section 7.9 we hand-coded a Newton iteration to solve the nonlinear BVP $$ 0.001u''-u^3 = 0,\qquad u(-1) = 1,~~ u(1) = -1. $$ Chebfun solves such problems automatically in response to the same syntax as above. Switching from L to N to suggest nonlinearity, let us write

N = chebop(@(x,u) 0.001*diff(u,2) - u^3);
N.lbc = 1; N.rbc = -1;

This gives us a chebop which Chebfun recognizes as nonlinear,

N =
   Nonlinear operator:
      u |--> 0.001.*diff(u,2)-u.^3
   operating on chebfun objects defined on:
    left boundary condition(s):
      u = 1
    right boundary condition(s):
      u = -1

To solve the BVP, we can write

u = N\0;
plot(u), grid on

Note that this is the same result as in Section 7.9. How does Chebfun solve such problems? That is a long story, which we shall not tell properly here. In brief, a Newton iteration (sometimes a damped Newton iteration) is carried out in "continuous mode", that is, in a space of functions rather than vectors. Recall that to find a zero of a scalar function, Newton's method requires a derivative at each iterative step, and to find a zero vector of a system of equations, it requires a Jacobian matrix. Here, we seek a zero function of a nonlinear differential operator equation. For this, Newton's method requires at each step the continuous analogue of a Jacobian matrix, which is a Frechet derivative linear operator. This Frechet derivative is realized in Chebfun by a continuous analogue of Automatic Differentiation using methods described in [Birkisson & Driscoll 2012].

Here is an example with a variable coefficient, a nonlinear BVP due to George Carrier analyzed in Sec. 9.7 of the book [Bender & Orzsag 1978]. We seek a function $u$ satisfying $$ \varepsilon u'' + 2(1-x^2) u + u^2 = 1, \qquad u(-1)=u(1) = 0 $$ with $\varepsilon = 0.01$. Here is a Chebfun formulation and solution.

ep = 0.01;
N = chebop(-1, 1);
N.op = @(x,u) ep*diff(u,2) + 2*(1 - x^2)*u + u^2;
N.bc = 'dirichlet';
u = N\1; plot(u), ylim([-2 2]), grid on

This is one of several valid solutions to this problem. To find another, we can specify a initial guess for the Newton iteration that differs from Chebfun's default (a polynomial function constructed to satisfy the boundary conditions --- the zero function in this case). For example, here we specify the initial guess $u(x) = 2(x^2 - 1)(1 - 2/(1 + 20x^2))$ and get a solution with two peaks instead of four.

x = chebfun('x');
N.init = 2*(x^2 - 1)*(1 - 2/(1 + 20*x^2));
[u, info] = solvebvp(N, 1);
plot(u), grid on

This time, instead of using \, we called the underlying method solvebvp, and we specified two output arguments. The second output is a MATLAB struct containing data showing the norms of the updates during the Newton iteration, revealing a slow initial phase followed by eventual rapid convergence.

nrmdu = info.normDelta;
semilogy(nrmdu,'.-k'), grid on, ylim([1e-14,1e2])

Another way to get information about the Newton iteration with nonlinear backlash is by setting




Type help cheboppref for details. Here we shall not pursue this option and thus return the system to its factory state:


When you apply backslash to a nonlinear chebop, it invokes the overloaded MATLAB command mldivide; this in turn calls solvebvp to do the actual work. By calling solvebvp directly, you can control the computation in ways not accessible through backslash, a situation just like the relationship between \ and linsolve for solving a linear system in MATLAB. See the help documentation for details.

10.2 Initial-value problems: \ and solveivp

For IVPs, Chebfun contains overloads ode113, ode45, and ode15s of familiar MATLAB codes. Again, however, these are not our recommended methods. Instead, we recommend \ and its equivalent solveivp.

For example, suppose we want to solve the nonlinear IVP $$ u' = u^2, \qquad t\in [0,1], \quad u(0)=0.95. $$ We can set up the problem like this:

N = chebop(0, 1);
N.op = @(t,u) diff(u) - u^2;
N.lbc = 0.95
N =
   Nonlinear operator:
      u |--> diff(u)-u.^2
   operating on chebfun objects defined on:
    left boundary condition(s):
      u = 0.95

Since boundary conditions have been specified at only one end of the domain, Chebfun knows that this is an initial value problem. We solve it and plot the solution:

u = N\0;
plot(u), grid on

A major change was introduced in Version 5.1 in how initial-value (and final-value) problems are solved. Before, Chebfun used the same global spectral representations as for BVPs. This usually works fine for linear problems, but for nonlinear ones, it is inferior to the method of time-stepping by Runge-Kutta or Adams formulas. In Chebfun Version 5.1, we accordingly switched to solving IVPs numerically by ode113 (by default), converting the resulting output to a Chebfun representation. (This work, a substantial job since higher-order equations must be reformulated as first-order systems, was carried out by Asgeir Birkisson [Brikisson 2018].) If you wish to invoke the global spectral method instead of time-stepping, you can write

u2 = solvebvp(N,0);

(With cheboppref.setDefaults('ivpSolver','collocation') you could make this switch globally.) For this problem the method converges, giving a solution that is close but not the same:

ans =

For many nonlinear IVPs, however, the solvebvp approach would not converge.

Here is an example of an IVP with two components; the equation happens to be linear.

N = chebop(0, 100);
N.op = @(t,u) diff(u,2) + u;
N.lbc = [1; 0];
u = N\0;
plot([u diff(u)]), grid on, ylim([-1.5 1.5])

As a third example, let us solve a van der Pol equation for a nonlinear oscillator, $$ \varepsilon u'' = (1-u^2)u' - u , \qquad t\in [0,20], ~~ u(0) = 3, ~~u'(0) = 0. $$ Here is a solution with $\varepsilon = 0.05$:

N = chebop(0,20);
N.op = @(t,u) 0.05*diff(u,2) - (1-u^2)*diff(u) + u;
N.lbc = [3; 0];
u = N\0;
plot(u), ylim([-4 4]), grid on

As a final example let us consider the famous Lorenz equations, whose solution trajectories are chaotic: We can set up the problem and solve it like this:

N = chebop(0,15);
N.op = @(t,u,v,w) [diff(u)-10*(v-u);
N.lbc = @(u,v,w) [u+14; v+15; w-20];
[u,v,w] = N\0;
plot3(u,v,w), view(-5,9), axis off

10.3 Stiff IVPs

Chebfun's default solution methods work well for moderately stiff ODE IVPs. For highly stiff problems, however, it is desirable to switch the underlying engine from the default ode113 to the stiff solver ode15s. For example, the problem $$ u' = -u - 10000(u(t)-\cos(t)), \qquad u(0) = 1 $$ has the solution $u(t) = \cos(t)$. However, it is highly stiff, and to solve it we can proceed as follows:

N = chebop(0, 10);
N.op = @(t,u) diff(u) + sin(t) + 10000*(u-cos(t));
N.lbc = 1;
pref = cheboppref; pref.ivpSolver = 'ode15s';
tic, u = solveivp(N,0,pref); toc
plot(u), grid on, ylim([-1.5 1.5])
Elapsed time is 0.329978 seconds.

If we don't specify ode15s, the solution takes minutes instead of seconds.

10.4 Periodic problems

A relatively new feature in Chebfun is the solution of periodic ODEs; see chapter 19 of [Trefethen, Birkisson & Driscoll 2018]. For example, here is a function encoded in the gallerytrig command:

plot(cheb.gallerytrig('tsunami'), 'color', [.6 .4 0]), grid on
ylim([-.2 .2])

If you look in gallerytrig, you will find that this curve has been generated by the following sequence:

op = @(x,u) diff(u,2) + diff(u) + 600*(1+sin(x))*u;
L = chebop(op, [-pi,pi], 'periodic');
f = L\1;

which corresponds to the problem $$ u'' + u' + 600(1+\sin(x)) = 0, \quad x\in [-\pi,\pi], ~~ u(-\pi)=u(\pi), ~~ u'(-\pi)=u'(\pi) . $$ The periodic flag instructs Chebfun to impose periodic boundary conditions and solve the problem with a trigonometric discretization, that is, a Fourier spectral method. The result is a trigfun:

f =
   chebfun column (1 smooth piece)
       interval       length     endpoint values trig
[    -3.1,     3.1]      109   -0.0043  -0.0043 
vertical scale = 0.17 

10.5 Ultraspherical discretizations

As with most Chebfun operations involving differential equations, for nonlinear ODE BVPs and periodic ODEs Chebfun offers a choice between the default spectral collocation methods or an alternative ultraspherical method. See Sections 7.7 and 8.10.

10.6 Graphical user interface: Chebgui

Chebfun includes a GUI (Graphical User Interface) called chebgui for interactive solution of ODE, time-dependent PDE, and eigenvalue problems. For many users, this is the single most important part of Chebfun. We will not describe chebgui here, but we encourage readers to give it a try. Be sure to note the Demo menu, which offers dozens of preloaded examples, both scalars and systems. Perhaps most important of all is the "Export to m-file" button, which produces a Chebfun m-file corresponding to whatever problem is loaded into the GUI. This feature enables one to get going quickly and interactively, then switch to a Chebfun program to adjust the fine points. To start exploring, just type chebgui.

10.7 References

[Bender & Orszag 1978] C. M. Bender and S. A. Orszag, Advanced Mathematical Methods for Scientists and Engineers, McGraw-Hill, 1978.

[Birkisson 2014] A. Birkisson, Numerical Solution of Nonlinear Boundary Value Problems for Ordinary Differential Equations in the Continuous Framework, D. Phil. thesis, University of Oxford, 2014.

[Birkisson 2018] A. Birkisson, Automatic reformulation of ODEs to systems of first order equations, Trans. Math. Softw., 44 (2018), 31.

[Birkisson & Driscoll 2012] A. Birkisson and T. A. Driscoll, Automatic Fréchet differentiation for the numerical solution of boundary-value problems, ACM Transactions on Mathematical Software, 38 (2012), 1-26.

[Birkisson & Driscoll 2013] A. Birkisson and T. A. Driscoll, Automatic linearity dection, preprint,, 2013.

[Trefethen, Birkisson & Driscoll 2013] L. N. Trefethen, A. Birkisson, and T. A. Driscoll, Exploring ODEs, SIAM, 2018; freely available at

% SIAM, 2018.