demo_matrixfunction2.m

% computes the solution at time t=1 of a heat equation
J = 30; h = 2/J; s = (-1+h:h:1-h)';          % 1D grid
[xx,yy] = meshgrid(s,s);                     % 2D grid
x = xx(:); y = yy(:);                        % 2D grid stretched to 1D
L = -gallery('poisson',J-1)/h^2;             % 2D Laplacian
u = (1-x.^2).*(1-y.^2).*exp(x);              % initial condition
v = cf_matrixexponential(0.02*L,u);          % exp(0.02*L)*u 

v_expm = expm(0.02*L)*u;
norm(v - v_expm)