function [x,y,t,u,Frames] = heat2d
% Heat diffusion in a slab

% Setup
L = 1.0;
W = 1.0;
T = 75;
K = 0.001;
n = 75;
dt = T/n;
l =10;
dx = L/l;
m = 10;
dy = W/m;

u = zeros(l+1,m+1,n+1);

% Create mesh
t = (0:n)*dt;
x = (0:l)*dx;
y = (0:m)*dy;

% Set i.c. and b.c.s
u(:,:,1) = u0(x,y)';

u(1,:,:) = bcLeft;

u(l+1,:,:) = bcRight;

u(:,1,:) = bcBottom;

u(:,m+1,:) = bcTop;




alphax = K*dt/dx^2;
alphay = K*dt/dy^2;

for k=1:n
    for j=2:m
        for i=2:l
            u(i,j,k+1) = alphax*u(i-1,j,k) + (1 - 2*alphax -2*alphay)*u(i,j,k)...
                + alphax*u(i+1,j,k) + alphay*u(i,j-1,k) + alphay*u(i,j+1,k);
        end
    end
end

 for k=1:n+1
     surf(x,y,u(:,:,k));
     axis([0 1 0 1 0 100]);
     view([0 90]);
     shading interp;
     Frames(k) = getframe;
 end

% Some functions

function u0 = u0(x,y)
u0 = 0.0;

function bcLeft = bcLeft
bcLeft = 100.0;

function bcRight = bcRight
bcRight = 0.0;

function bcBottom = bcBottom
bcBottom = 100.0;

function bcTop = bcTop
bcTop = 0.0;
