% 2D Heat Equation using discretization:
clear all; 
L=20; N=64; tspan=[0:0.1:10];
x2=linspace(-L/2,L/2,N+1); x=x2(1:N); dx=x(2)-x(1);
y2=linspace(-L/2,L/2,N+1); y=y2(1:N);
[X,Y]=meshgrid(x,y);
%
kappa = .1;
u0=exp(-(X.^2+Y.^2));

e1=ones(N,1);
A=spdiags([e1 -2*e1 e1],[-1 0 1],N,N);
A=A/(dx^2);
I=eye(N);
B=kron(I,A)+kron(A,I);
U0 = reshape(u0,N^2,1);
[t,y]=ode45('heat_2d_rhs',tspan,U0,[],B,N,kappa);
M = moviein(length(tspan));
for j=1:length(tspan)
  u=reshape(y(j,1:N^2),N,N);
  surf(X,Y,u), shading interp
  axis([-L/2 L/2 -L/2 L/2 0 1]);
  xlabel('x'); ylabel('y'); zlabel('u');
  M(:,j) = getframe;
end