% Simple Heat Equation Solver
% 5-21-2008

clear all; close all; clc

% Problem setup
sigma = 1; k = 1;
N=100;
L=1;
tspan = [0 0.4];
x=linspace(-L,L,N);
dx=mean(diff(x));

% Construct the matrix A
e1=ones(N,1);
A = spdiags([e1 -2*e1 e1],[-1 0 1], N,N);
A(1,N) = 1;  A(N,1) = 1;

% Generate u0
u0 = 1/(sigma*sqrt(2*pi)) * exp(-x.^2/sigma^2);

[t,u] = ode45('heat_rhs',tspan,u0,[],k,dx,A);

surf(x,t,u); xlabel('x'); ylabel('t'); zlabel('u');
shading flat;