% heatmexp.m
% solves heat equation 
% using matrix exponential

% Example:

% [x,y,y0,ya] = heatmexp(121);
% plot(x,y0,x,y,x,ya,'o')

function [x,y,y0,ya] = heatmexp(n)

% initial condition

x = linspace(-10,10,n)';

dxs = (x(2)-x(1))^2;

y0 = 1/sqrt(2*pi)*exp(-x.^2/2);

y = y0;

tMax = 1;
nSteps = 1000;
dt = tMax/nSteps;

% set up linear propagator: matrix exponential

M = -2*diag(ones(1,n)) + diag(ones(1,n-1),1) + diag(ones(1,n-1),-1);

kappa = 0.5;

L = expm(M/dxs*kappa*dt);

for k=1:nSteps
	
    y = L*y;
    
end

% done with numerics, compute
% analytical solution

ya = 1/sqrt(2*pi*2)*exp(-x.^2/4);

end