1.1 What is a chebfun?

A chebfun is a function of one variable defined on an interval $[a,b]$. The syntax for chebfuns is almost exactly the same as the usual MATLAB syntax for vectors, with the familiar MATLAB commands for vectors overloaded in natural ways. Thus, for example, whereas sum(f) returns the sum of the entries when f is a vector, it returns a definite integral when f is a chebfun.

Chebfun with a capital C is the name of the software system.

The aim of Chebfun is to "feel symbolic but run at the speed of numerics". More precisely, our vision is to achieve for functions what floating-point arithmetic achieves for numbers: rapid computation in which each successive operation is carried out exactly apart from a rounding error that is very small in relative terms [Trefethen 2007].

The implementation of Chebfun is based on the mathematical fact that smooth functions can be represented very efficiently by polynomial interpolation in Chebyshev points, or equivalently, thanks to the Fast Fourier Transform, by expansions in Chebyshev polynomials. For a simple function, 20 or 30 points often suffice, but the process is stable and effective even for functions complicated enough to require 1000 or 1,000,000 points. Chebfun makes use of adaptive procedures that aim to find the right number of points automatically so as to represent each function to roughly machine precision, that is, about 15 digits of relative accuracy. (Originally Chebfun stored function values at Chebyshev points; in Version 5 it switched to storing Chebyshev expansion coefficients.)

The mathematical foundations of Chebfun are for the most part well established by results scattered throughout the 20th century. A key early figure, for example, was Bernstein in the 1910s. Much of the relevant material can be found collected in the Chebfun-based book Approximation Theory and Approximation Practice [Trefethen 2013].

Chebfun was originally created by Zachary Battles and Nick Trefethen at Oxford during 2002-2005 [Battles & Trefethen 2004]. Battles left the project in 2005, and soon four new members were added to the team: Ricardo Pachon (from 2006), Rodrigo Platte (from 2007), and Toby Driscoll and Nick Hale (from 2008). In 2009, Asgeir Birkisson and Mark Richardson also became involved, and other contributors included Pedro Gonnet, Joris Van Deun, and Georges Klein. Nick Hale served as Director of the project during 2010-2014. The Chebfun Version 5 rewrite was directed by Nick Hale during 2013-2014, and the team included Anthony Austin, Asgeir Birkisson, Toby Driscoll, Hrothgar, Mohsin Javed, Hadrien Montanelli, Alex Townsend, Nick Trefethen, Grady Wright, and Kuan Xu. October 2014 brough new arrivals Jared Aurentz, Behnam Hashemi, and Mikael Slevinsky. Further information about Chebfun history is available at the Chebfun web site, http://www.chebfun.org where one can also find a discussion of other software projects related to Chebfun. This Guide is based on Chebfun Version 5.3, released in December 2015.

1.2 Constructing simple chebfuns

The chebfun command constructs a chebfun from a specification such as a string or an anonymous function. If you don't specify an interval, then the default interval $[-1,1]$ is used. For example, the following command makes a chebfun corresponding to $\cos(20x)$ on $[-1,1]$ and plots it.

  f = chebfun('cos(20*x)');
  plot(f), ylim([-1.2,1.2])

From this little experiment, you cannot see that f is represented by a polynomial. One way to see this is to find the length of f:

  length(f)
ans =
    51

Another is to remove the semicolon that suppresses output:

  f
f =
   chebfun column (1 smooth piece)
       interval       length     endpoint values  
[      -1,       1]       51      0.41     0.41 
vscale = 1.000000e+00.

These results tell us that f is represented by a polynomial interpolant through 51 Chebyshev points, i.e., a polynomial of degree 50. These numbers have been determined by an adaptive process. We can see the data points by plotting f with the '.-' option:

  plot(f,'.-'), ylim([-1.2 1.2])

The formula for $N+1$ Chebyshev points in $[-1,1]$ is $$ x(j) = -\cos(j \pi/N), \quad j = 0:N, $$ and in the figure we can see that the points are clustered accordingly near $1$ and $-1$. Note that in the middle of the grid, there are about 5 points per wavelength, which is evidently what it takes to represent this cosine to 15 digits of accuracy. For intervals other than $[-1,1]$, appropriate Chebyshev points are obtained by a linear scaling.

The curve between the data points is the polynomial interpolant, which can be evaluated by the barycentric formula introduced by Salzer [Berrut & Trefethen 2004, Salzer 1972]. This method of evaluating polynomial interpolants is stable and efficient even if the degree is in the millions [Higham 2004]. In recent years Chebfun actually evaluates polynomials from their Chebyshev series rather than by barycentric interpolation; the difference in the two methods is little.

What is the integral of $f$ from $-1$ to $1$? Here it is:

  sum(f)
ans =
   0.091294525072762

This number was computed by integrating the polynomial (Clenshaw-Curtis quadrature -- see Section 2.1), and it is interesting to compare it to the exact answer from calculus:

  exact = sin(20)/10
exact =
   0.091294525072763

Here is another example, now with the chebfun defined by an anonymous function instead of a string. In this case the interval is specified as $[0,100]$.

  g = chebfun(@(t) besselj(0,t),[0,100]);
  plot(g), ylim([-.5 1])

The function looks complicated, but it is actually a polynomial of surprisingly small degree:

  length(g)
ans =
    89

Is it accurate? Well, here are three random points in $[0,100]$:

  format long
  x = 100*rand(3,1)
x =
  81.472368639317892
  90.579193707561927
  12.698681629350606

Let's compare the chebfun to the true Bessel function at these points:

  exact = besselj(0,x);
  error = g(x) - exact;
  [g(x) exact error]
ans =
   0.048059538377880   0.048059538377880   0.000000000000000
  -0.021311086924089  -0.021311086924090   0.000000000000000
   0.176415464952875   0.176415464952875  -0.000000000000000

If you want to know the first 5 zeros of the Bessel function, here they are:

  r = roots(g); r = r(1:5)
r =
   2.404825557695767
   5.520078110286308
   8.653727912911007
  11.791534439014278
  14.930917708487790

Notice that we have just done something nontrivial and potentially useful. How else would you find zeros of the Bessel function so readily? As always with numerical computation, we cannot expect the answers to be exactly correct, but they will usually be very close. In fact, these computed zeros are accurate to close to machine precision:

  besselj(0,r)
ans =
   1.0e-14 *
   0.314549201444184
  -0.082399534411417
   0.154471710674566
  -0.058908346770306
  -0.092296985262750

Most often we get a chebfun by operating on other chebfuns. For example, here is a sequence that uses plus, times, divide, and power operations on an initial chebfun x to produce a famous function of Runge:

  x = chebfun('x');
  f = 1./(1+25*x.^2);
  length(f)
  clf, plot(f)
ans =
   181

1.3 Operations on chebfuns

There are more than 200 commands that can be applied to a chebfun. For a list of many of them you can type methods:

  methods chebfun
Methods for class chebfun:

abs              csc              jump             rank             
acos             cscd             kron             rdivide          
acosd            csch             ldivide          real             
acosh            ctranspose       le               reallog          
acot             cumprod          legcoeffs        realpow          
acotd            cumsum           legpoly          realsqrt         
acoth            cylinder         length           rem              
acsc             diff             log              remez            
acscd            dirac            log10            removeDeltas     
acsch            disp             log1p            repmat           
airy             display          log2             residue          
all              domain           logical          restrict         
and              eigs             loglog           roots            
angle            ellipj           lt               round            
any              ellipke          lu               sec              
arcLength        end              mat2cell         secd             
area             eq               max              sech             
asec             erf              mean             semilogx         
asecd            erfc             measure          semilogy         
asech            erfcinv          merge            sign             
asin             erfcx            mesh             simplify         
asind            erfinv           min              sin              
asinh            exp              minandmax        sinc             
atan             expm             minus            sind             
atan2            expm1            mldivide         sinh             
atan2d           feval            mod              size             
atand            fill             movie            sound            
atanh            find             mrdivide         spy              
besselh          fix              mtimes           sqrt             
besseli          fliplr           nchoosek         std              
besselj          flipud           ne               subsasgn         
besselk          floor            newDomain        subspace         
bessely          fracDiff         nextpow2         subsref          
bvp4c            fracInt          norm             sum              
bvp5c            fred             normal           surf             
cat              gamma            normest          surface          
ceil             ge               not              surfc            
cf               get              null             svd              
cheb2cell        gmres            num2cell         tan              
cheb2quasi       gt               or               tand             
chebcoeffs       heaviside        orth             tanh             
chebellipseplot  horzcat          overlap          times            
chebfun          hscale           pde15s           transpose        
chebpade         hypot            pde23t           trigcoeffs       
chebpoly         imag             permute          trigremez        
chebtune         innerProduct     pinv             truncate         
circconv         integral         plot             uminus           
comet            inv              plot3            unwrap           
comet3           isPeriodicTech   plotcoeffs       uplus            
complex          isdelta          plotregion       vander           
compose          isempty          plus             var              
cond             isequal          poly             vertcat          
conj             isfinite         polyfit          volt             
conv             ishappy          polyval          vscale           
cos              isinf            pow2             waterfall        
cosd             isnan            power            why              
cosh             isreal           prod             xor              
cot              issing           qr               
cotd             iszero           quantumstates    
coth             jaccoeffs        quasi2cheb       
cov              join             range            

Static methods:

dct              idlt             ndct             pchip            
dlt              idst             ode113           spline           
dst              interp1          ode15s           update           
idct             lagrange         ode45            

To find out what a command does, you can use help.

  help chebfun/sum
 SUM   Definite integral of a CHEBFUN.
    SUM(F) is the integral of a column CHEBFUN F over its domain of definition.

    SUM(F, A, B), where A and B are scalars, integrates a column CHEBFUN F over
    [A, B], which must be a subdomain of F.domain:

                          B
                          /
                SUM(F) =  | F(t) dt.
                          /
                         A

    SUM(F, A, B), where A and B are CHEBFUN objects, returns a CHEBFUN S which
    satisfies

                        B(s)
                        /
                S(s) =  | F(t) dt.
                        /
                      A(s)

    SUM(F, DIM), where DIM is one of 1, 2, sums F over the dimension DIM. If F
    is a column CHEBFUN and DIM = 1 or if F is a row CHEBFUN and DIM = 2 then
    this integrates in the continuous dimension of F, as described above.
    Otherwise, SUM(F, DIM) sums across the columns (rows) of the column (row)
    CHEBFUN F.

  See also CUMSUM, DIFF.

Most of the commands in the list above exist in ordinary MATLAB; some exceptions are domain, restrict, chebpoly, and remez. We have already seen length and sum in action. In fact we have already seen subsref too, since that is the MATLAB command for (among other things) evaluating arguments in parentheses. Here is another example of its use:

  f(0.5)
ans =
   0.137931034482759

Here for comparison is the true result:

  1/(1+25/4)
ans =
   0.137931034482759

In this Runge function example, we have also implicitly seen times, plus, power, and rdivide, all of which have been overloaded from their usual MATLAB uses to apply to chebfuns.

In the next part of this tour we shall explore many of these commands systematically. First, however, we should see that chebfuns are not restricted to smooth functions.

1.4 Piecewise smooth chebfuns

Many functions of interest are not smooth but piecewise smooth. In this case a chebfun may consist of a concatenation of smooth pieces, each with its own polynomial representation. Each of the smooth pieces is called a "fun". This enhancement of Chebfun was developed initially by Ricardo Pachon during 2006-2007, then also by Rodrigo Platte starting in 2007 [Pachon, Platte and Trefethen 2010]. Essentially funs are the "classic chebfuns" for smooth functions on $[-1,1]$ originally implemented by Zachary Battles in Chebfun Version 1.

Later we shall describe the options in greater detail, but for the moment let us see some examples. One way to get a piecewise smooth function is directly from the Chebfun constructor, taking advantage of its capability of automatic edge detection. For example, in the default "splitting off" mode a function with a jump in its derivative produces a warning message,

  f = chebfun('abs(x-.3)');
Warning: Function not resolved using 65537 pts. Have you tried 'splitting on'? 

The same function can be successfully captured with splitting on:

  f = chebfun('abs(x-.3)','splitting','on');

The length command reveals that f is defined by four data points, two for each linear interval:

  length(f)
ans =
     4

We can see the structure of f in more detail by typing f without a semicolon:

  f
f =
   chebfun column (2 smooth pieces)
       interval       length     endpoint values  
[      -1,     0.3]        2       1.3  3.3e-16 
[     0.3,       1]        2   1.1e-16      0.7 
vscale = 1.300000e+00.   Total length = 4

This output confirms that f consists of two funs, each defined by two points and two corresponding function values. The functions live on intervals defined by breakpoints at $-1$, $1$, and a number very close to $0.3$.

Another way to make a piecewise smooth chebfun is to construct it explicitly from various pieces. For example, the following command specifies three functions $x^2$, $1$, and $4-x$, together with a vector of endpoints indicating that the first function applies on $[-1,1]$, the second on $[1,2]$, and the third on $[2,4]$:

  f = chebfun({@(x) x.^2, 1, @(x) 4-x},[-1 1 2 4]);
  plot(f)

We expect f to consist of three pieces of lengths 3, 1, and 2, and this is indeed the case:

  f
f =
   chebfun column (3 smooth pieces)
       interval       length     endpoint values  
[      -1,       1]        3         1        1 
[       1,       2]        1         1        1 
[       2,       4]        2         2 -1.1e-16 
vscale = 2.   Total length = 6

Our eyes see pieces, but to Chebfun, f is just another function. For example, here is its integral.

  sum(f)
ans =
   3.666666666666667

Here is an algebraic transformation of f, which we plot in another color for variety.

  plot(1./(1+f),'r')

Some Chebfun commands naturally introduce breakpoints in a chebfun. For example, the abs command first finds zeros of a function and introduces breakpoints there. Here is a chebfun consisting of 6 funs:

  f = abs(exp(x).*sin(8*x));
  plot(f)

And here is an example where breakpoints are introduced by the max command, leading to a chebfun with 13 pieces:

  f = sin(20*x);
  g = exp(x-1);
  h = max(f,g);
  plot(h), ylim([0 1.2])

As always, h may look complicated to a human, but to Chebfun it is just a function. Here are its mean, standard deviation, minimum, and maximum:

  mean(h)
ans =
   0.578242020778010
  std(h)
ans =
   0.280937455806246
  min(h)
ans =
   0.135335283236613
  max(h)
ans =
   1.000000000000000

A final note about piecewise smooth chebfuns is that the automatic edge detection or "splitting" feature, when it is turned on, may subdivide functions even though they do not have clean point singularities, and this may be desirable or undesirable depending on the application. For example, considering $\sin(x)$ over $[0,1000\pi]$ with splitting on, we end up with a chebfun with many pieces:

  tic, f = chebfun('sin(x)',[0 1000*pi],'splitting','on'), toc
f =
   chebfun column (32 smooth pieces)
       interval       length     endpoint values  
[       0,      98]       91   3.3e-14    -0.71 
[      98,   2e+02]       88     -0.71        1 
[   2e+02, 2.9e+02]       88         1    -0.71 
[ 2.9e+02, 3.9e+02]       87     -0.71  1.7e-13 
[ 3.9e+02, 4.9e+02]       87   3.1e-13     0.71 
[ 4.9e+02, 5.9e+02]       87      0.71       -1 
[ 5.9e+02, 6.9e+02]       86        -1     0.71 
[ 6.9e+02, 7.9e+02]       85      0.71 -7.7e-13 
[ 7.9e+02, 8.8e+02]       86  -1.3e-12    -0.71 
[ 8.8e+02, 9.8e+02]       86     -0.71        1 
[ 9.8e+02, 1.1e+03]       86         1    -0.71 
[ 1.1e+03, 1.2e+03]       86     -0.71 -3.1e-13 
[ 1.2e+03, 1.3e+03]       85  -6.6e-13     0.71 
[ 1.3e+03, 1.4e+03]       86      0.71       -1 
[ 1.4e+03, 1.5e+03]       86        -1     0.71 
[ 1.5e+03, 1.6e+03]       85      0.71 -1.6e-12 
[ 1.6e+03, 1.7e+03]       86  -1.3e-12    -0.71 
[ 1.7e+03, 1.8e+03]       86     -0.71        1 
[ 1.8e+03, 1.9e+03]       86         1    -0.71 
[ 1.9e+03,   2e+03]       85     -0.71  8.3e-13 
[   2e+03, 2.1e+03]       85   2.2e-12     0.71 
[ 2.1e+03, 2.2e+03]       86      0.71       -1 
[ 2.2e+03, 2.3e+03]       86        -1     0.71 
[ 2.3e+03, 2.4e+03]       87      0.71 -1.7e-12 
[ 2.4e+03, 2.5e+03]       85   1.5e-12    -0.71 
[ 2.5e+03, 2.6e+03]       86     -0.71        1 
[ 2.6e+03, 2.7e+03]       84         1    -0.71 
[ 2.7e+03, 2.7e+03]       85     -0.71 -3.2e-12 
[ 2.7e+03, 2.8e+03]       85  -2.5e-13     0.71 
[ 2.8e+03, 2.9e+03]       84      0.71       -1 
[ 2.9e+03,   3e+03]       86        -1     0.71 
[   3e+03, 3.1e+03]       85      0.71 -3.7e-12 
vscale = 1.000000e+00.   Total length = 2752
Elapsed time is 0.621162 seconds.

In this case it is more efficient -- and more interesting mathematically -- to omit the splitting and construct one global chebfun:

  tic, f2 = chebfun('sin(x)',[0 1000*pi]), toc
f2 =
   chebfun column (1 smooth piece)
       interval       length     endpoint values  
[       0, 3.1e+03]     1684  -3.4e-16 -3.2e-13 
vscale = 9.999862e-01.
Elapsed time is 0.040536 seconds.

Splitting on and off are discussed further in Section 8.3.

1.5 Infinite intervals and infinite function values

A major change from Chebfun Version 2 to Version 3 was the generalization of chebfuns to allow certain functions on infinite intervals or which diverge to infinity; the initial credit for these innovations belongs to Nick Hale, Rodrigo Platte, and Mark Richardson, and many later developments are due to Kuan Xu. For example, here is a function on the whole real axis,

  f = chebfun('exp(-x.^2/16).*(1+.2*cos(10*x))',[-inf,inf]);
  plot(f)

and here is its integral:

  sum(f)
Warning: Result may not be accurate as the function decays slowly at infinity. 
ans =
   7.089815392600106

Here's the integral of a function on $[1,\infty)$:

  sum(chebfun('1./x.^4',[1 inf]))
ans =
   0.333333333328450

Notice that several digits of accuracy have been lost here. Be careful! -- operations involving infinities in Chebfun are not always as accurate and robust as their finite counterparts.

Here is an example of a function that diverges to infinity, which we can capture with the 'exps' flag; see Chapter 7 for details:

  h = chebfun('(1/pi)./sqrt(1-x.^2)','exps',[-.5 -.5]);
  plot(h)

In this case the integral comes out just right:

  sum(h)
ans =
   1.000000000000000

For more on the treatment of infinities in Chebfun, see Chapter 9.

1.6 Periodic functions

Until 2014, Chebfun used only nonperiodic representations, based on Chebyshev polynomials. Beginning with Version 5, there is a new capability of representing sufficiently smooth periodic functions by trigonometric polynomials instead, that is, Fourier series. Such an object is still called a chebfun, but it is a periodic one, and the signal to invoke such capabilities is the string trig. For abbreviation, we call a periodic chebfun a ``trigfun''. This section gives a quick introduction, and more details can be found in Chapter 11.

Trigfuns were initiated by Grady Wright in the first half of 2014 [Wright et al. 2015]. A very interesting project along the same lines was carried out independently by Kristyn McLeod and Rodrigo Platte at Arizona State University [McLeod 2014].

For example, here is a periodic function on $[-\pi,\pi]$ represented in the usual way by a Chebyshev series.

ff = @(t) sin(t) + cos(2*t) - cos(t)/3 + cos(100*t)/6;
f = chebfun(ff,[-pi,pi]);
max(f)
plot(f)
ans =
   1.614526099978749

Its length, very roughly, is $100 \pi$,

length(f)
ans =
   383

Here is the same function represented by a Fourier series:

f2 = chebfun(ff,[-pi,pi],'trig')
max(f2)
plot(f2,'m')
f2 =
   chebfun column (1 smooth piece)
       interval       length     endpoint values trig
[    -3.1,     3.1]      201       1.5      1.5 
vscale = 2.122902e+00.
ans =
   1.614526099978751

Its length is now only about $200$ (exactly 201). This improvement by a factor of about $\pi/2$ is typical.

length(f2)
ans =
   201

Sampling at a few arbitrary points confirms that the two functions agree closely:

xx = [1/3 sqrt(2) exp(1)];
f(xx) - f2(xx)
ans =
   1.0e-14 *
   0.510702591327572   0.954791801177635   0.710542735760100

Readers may be interested to compare plotcoeffs applied to the first and second versions of $f$. Rather than display that here we shall turn to a simpler example involving a shorter Fourier series. Consider the function

f = chebfun('7 + sin(t) + exp(1)*cos(2*t)',[-pi,pi],'trig')
f =
   chebfun column (1 smooth piece)
       interval       length     endpoint values trig
[    -3.1,     3.1]        5       9.7      9.7 
vscale = 9.718282e+00.

Here are the coefficients of $f$ as an expansion in sines and cosines:

[a,b] = trigcoeffs(f)
a =
   7.000000000000000
   0.000000000000000
   2.718281828459045
b =
   1.000000000000000
   0.000000000000000

Here they are as an expansion in complex exponentials:

c = trigcoeffs(f)
c =
  1.359140914229523 + 0.000000000000000i
  0.000000000000000 + 0.500000000000000i
  7.000000000000000 + 0.000000000000000i
  0.000000000000000 - 0.500000000000000i
  1.359140914229523 - 0.000000000000000i

Bookkeeping of Fourier coefficients can often be a headache. If these examples don't make the patterns clear, details can be found with help trigcoeffs.

For a mathematically less trivial example, here is the cosine expansion of a function whose Fourier series coefficients are known to be values of a Bessel function:

f = chebfun('exp(cos(t))',[-pi pi],'trig');
[a,b] = trigcoeffs(f);
n = floor(length(f)/2);
exact = 2*besseli(0:n,1); exact(1) = exact(1)/2;
disp('        computed             exact')
disp([a exact'])
        computed             exact
   1.266065877752008   1.266065877752008
   1.130318207984970   1.130318207984970
   0.271495339534077   0.271495339534077
   0.044336849848664   0.044336849848664
   0.005474240442094   0.005474240442094
   0.000542926311914   0.000542926311914
   0.000044977322954   0.000044977322954
   0.000003198436462   0.000003198436462
   0.000000199212481   0.000000199212481
   0.000000011036772   0.000000011036772
   0.000000000550590   0.000000000550590
   0.000000000024980   0.000000000024980
   0.000000000001039   0.000000000001039
   0.000000000000040   0.000000000000040
   0.000000000000001   0.000000000000001

1.7 Rows, columns, and quasimatrices

MATLAB doesn't only deal with column vectors: there are also row vectors and matrices. The same is true of Chebfun. The chebfuns shown so far have all been in column orientation, which is the default, but one can also take the transpose, compute inner products, and so on:

  x = chebfun(@(x) x)
x =
   chebfun column (1 smooth piece)
       interval       length     endpoint values  
[      -1,       1]        2        -1        1 
vscale = 1.
  x'
ans =
   chebfun row (1 smooth piece)
       interval       length     endpoint values  
[      -1,       1]        2        -1        1 
vscale = 1.
  x'*x
ans =
   0.666666666666667

One can also make matrices whose columns are chebfuns or whose rows are chebfuns, like this:

  A = [1 x x.^2]
A =
   chebfun column1 (1 smooth piece)
       interval       length     endpoint values  
[      -1,       1]        3         1        1 
vscale = 1.
   chebfun column2 (1 smooth piece)
       interval       length     endpoint values  
[      -1,       1]        3        -1        1 
vscale = 1.
   chebfun column3 (1 smooth piece)
       interval       length     endpoint values  
[      -1,       1]        3         1        1 
vscale = 1.
  A'*A
ans =
   2.000000000000000  -0.000000000000000   0.666666666666667
  -0.000000000000000   0.666666666666667   0.000000000000000
   0.666666666666667   0.000000000000000   0.400000000000000

These are called quasimatrices, and they are discussed in Chapter 6.

1.8 Chebfun features not in this Guide

Some of Chebfun's most remarkable features haven't made it into this edition of the Guide. Here are some of our favorites:

o leg2cheb and cheb2leg for fast Legendre-Chebyshev conversions (also legvals2chebcoeffs, chebcoeffs2legpts, and ten more)

o conv and circconv for convolution,

o The 'equi' flag to the Chebfun constructor for equispaced data,

o polyfit for least-squares fitting in the continuous context,

o inv for computing the inverse of a chebfun,

o pde15s for PDEs in one space and one time variable.

To learn about any of these options, try the appropriate help command. Just as a taster, here's a hint of how fast Chebfun can convert a ten-thousand coefficient Chebyshev expansion to Legendre coefficients and back again using an algorithm from [Hale & Townsend 2013]:

tic
ccheb = randn(10000,1);
cleg = cheb2leg(ccheb);
ccheb2 = leg2cheb(cleg);
norm(ccheb-ccheb2,inf)
toc
ans =
     3.814496635223819e-11
Elapsed time is 0.457952 seconds.

1.9 Chebfun example galleries

MATLAB has long had a gallery command to generate interesting matrices, and with version 5.1, Chebfun has introduced an analogous gallery command to generate interesting functions.

Here is what is currently available:

help cheb.gallery
 GALLERY   Chebfun example functions.
    F = GALLERY(NAME) returns a chebfun or a quasimatrix corresponding to NAME.
    See the listing below for available names.

    For example,  plot(cheb.gallery('zigzag'))  plots a degree 10000 polynomial
    that doesn't look like a polynomial, and  plot(cheb.gallery('gamma'))  shows
    a chebfun with poles. For details of how each function is constructed, try
    type +cheb/gallery  or  edit cheb.gallery.

    [F,FA] = GALLERY(NAME) also returns the anonymous function FA used to
    define the function. Some gallery functions are generated by operations
    beyond the usual Chebfun constructor (e.g. by solving ODEs), so FA in those
    cases simply evaluates the chebfun.

    GALLERY with no input argument returns a random function from the gallery.

    GALLERY with no output argument creates a plot of the selected function.

    airy         Airy Ai function on [-40,40]
    bessel       Bessel function J_0 on [-100,100]
    bump         C-infinity function with compact support
    blasius      Blasius function on [0,10]
    chirp        Sine with exponentially increasing frequency
    daubechies   Approximation to Daubechies phi_2 wavelet scaling function
    erf          Error function on [-10,10]
    fishfillet   Wild oscillations from Extreme Extrema example
    gamma        Gamma function on [-4,4]
    gaussian     Gaussian function on [-Inf,Inf]
    jitter       A piecewise constant function generated by ROUND
    kahaner      Challenging integrand with four spikes
    motto        Chebfun motto (Gilbert Strang)
    random       Polynomial interpolant through random data in Chebyshev points
    rose         A complex-valued sinusoid
    runge        Runge function
    seismograph  Tanh plus growing oscillation
    Si           Sine integral on [-50,50]
    sinefun1     As smooth as it looks
    sinefun2     Not as smooth as it looks
    spikycomb    25 peaks, each sharper than the last
    stegosaurus  max(wiggly, x/10)
    vandercheb   Chebyshev-Vandermonde quasimatrix
    vandermonde  Vandermonde quasimatrix
    wiggly       One of the Chebfun team's favorites
    zigzag       Degree 10000 polynomial that looks piecewise linear

    Gallery functions are subject to change in future releases of Chebfun.

  See also CHEB.GALLERYTRIG, CHEB.GALLERY2.

For example, here is a chebfun representing the Airy function,

plot(cheb.gallery('airy')), ylim([-.8 .8])
title('Airy function')

In this instance the underlying code fits in a line,

fa = @airy; p = chebfun(fa, [-40,40]);

Some examples make use of more complicated code, like this approximation to a Daubechies wavelet scaling function (accurate to about 3 digits of accuracy; the underlying function is a fractal):

plot(cheb.gallery('daubechies')), ylim([-0.5 1.5])
title('Daubechies scaling function')

To find out how a gallery example was generated, take a look at the code with type +cheb/gallery or edit +cheb/gallery.

Like the MATLAB gallery command, cheb.gallery produces a plot if you call it without specifying output variables. To illustrate, let us finish with an example the Chebfun team enjoys from the appendix to [Trefethen 2013], "Six myths of polynomial interpolation and quadrature":

cheb.gallery('zigzag')

This function looks piecewise linear, but in fact, it is a polynomial of degree 10000. This serves no purpose from an approximation point of view -- one would never represent this function in this manner -- but it illustrates the robustness of high-degree polynomial approximation.

If you call cheb.gallery without any input arguments, it selects a gallery function at random.

Other collections worth exploring are cheb.gallerytrig for periodic functions and cheb.gallery2 for 2D functions.

1.10 How this Guide is produced

This guide is produced in MATLAB using the publish command with a style sheet somewhat different from the usual; the output of publish is then processed by Markdown. To publish a chapter for yourself, make sure the chebfun guide directory is in your path and then type, for example, open(publish('guide1')). The formatting may not be exactly right but it should certainly be intelligible.

1.11 References

[Battles & Trefethen 2004] Z. Battles and L. N. Trefethen, "An extension of MATLAB to continuous functions and operators", SIAM Journal on Scientific Computing, 25 (2004), 1743-1770.

[Berrut & Trefethen 2005] J.-P. Berrut and L. N. Trefethen, "Barycentric Lagrange interpolation", SIAM Review 46, (2004), 501-517.

[Hale & Townsend 2013] N. Hale and A. Townsend, A fast, simple, and stable Chebyshev--Legendre transform using an asymptotic formula, SIAM Journal on Scientific Computing, 36 (2014), A148-A167.

[Higham 2004] N. J. Higham, "The numerical stability of barycentric Lagrange interpolation", IMA Journal of Numerical Analysis, 24 (2004), 547-556.

[McLeod 2014] K. N. McLeod, "Fourfun: A new system for automatic computations using Fourier expansions," SIAM Undergraduate Research Online, 7 (2014), http://dx.doi.org/10.1137/14S013238.

[Pachon, Platte & Trefethen 2010] R. Pachon, R. B. Platte and L. N. Trefethen, "Piecewise-smooth chebfuns", IMA J. Numer. Anal., 30 (2010), 898-916.

[Salzer 1972] H. E. Salzer, "Lagrangian interpolation at the Chebyshev points cos(nu pi/n), nu = 0(1)n; some unnoted advantages", Computer Journal 15 (1972), 156-159.

[Trefethen 2007] L. N. Trefethen, "Computing numerically with functions instead of numbers", Mathematics in Computer Science 1 (2007), 9-19. Revised and reprinted in Communications of the ACM 58 (2014), 91-97.

[Trefethen 2013] L. N. Trefethen, Approximation Theory and Approximation Practice, SIAM, 2013.

[Wright et al. 2015] G. B. Wright, M. Javed, H. Montanelli, and L. N. Trefethen, Extension of Chebfun to periodic functions, SIAM J. Sci. Comp., to appear.