A well-known PDE problem is the heat equation initial boundary-value problem posed for $x\in [a, b]$ and $t>0$,

$$ u_t = u_{xx}, ~~~~ u(x, 0) = u_0(x) $$

with suitable boundary conditions. We can regard this as a time-dependent linear process

$$ u_t = Lu $$

where $L$ is the operator $d^2/dx^2$ on $[a, b]$ with the same boundary conditions. The solution is

$$ u(t) = \exp(tL) u(0). $$

In Chebfun we can implement this idea using the expm command to compute the operator exponential. Here is an example with Neumann boundary conditions on the interval $[0, 6]$. We start with quite an irregular initial function.

d = [0, 6];
u = chebfun(@(x) sign((-1).^floor(x.^1.5)), d, 'splitting', 'on');
LW = 'linewidth'; lw = 2; FS = 'fontsize'; fs = 16;
clf, plot(u, LW, lw), grid on
title(sprintf('t = %4.2f     length = %d', 0, length(u)), FS, fs)
ax = [0 6 -1.2 1.2]; axis(ax)

Here's the solution at $t = 0.01$. Notice that the narrower spikes have lost more amplitude than the wider ones.

L = chebop(d);            % operator on domain [0, 6]
L.op = @(u) diff(u, 2);   % 2nd-derivative operator
L.lbc = @(u) diff(u);     % Neumann BC at left
L.rbc = @(u) diff(u);     % Neumann BC at right
dt = 0.01;
u = expm(L, dt, u);       % exponential of the operator
plot(u, LW, lw), axis(ax), grid on
title(sprintf('t = %4.2f     length = %d', 0.01, length(u)), FS, fs)

Here is the solution at $t = 0.02$. The rightmost maximum has extra amplitude, since it effectively corresponded to a wider initial spike thanks to the Neumann boundary condition.

u = expm(L, dt, u);
plot(u, LW, lw), axis(ax), grid on
title(sprintf('t = %4.2f     length = %d', 0.02, length(u)), FS, fs)

At $t=0.1$, there is not much of the original structure left. The length of the chebfun has also been reduced.

u = expm(L, 8*dt, u);
plot(u, LW, lw), axis(ax), grid on
title(sprintf('t = %4.2f     length = %d', 0.1, length(u)), FS, fs)