2d heat equation in matlab

main.m
% Define the domain and initial conditions
Lx = 1; % length of x-domain
Ly = 1; % length of y-domain
x = linspace(0, Lx, 100);
y = linspace(0, Ly, 100);
[x, y] = meshgrid(x, y);
u0 = sin(pi*x).*sin(pi*y);

% Define the diffusion coefficient and the boundary conditions
D = 1; % diffusion coefficient
xmesh = linspace(0,Lx,10);
ymesh = linspace(0,Ly,10);
msh = createMesh2D(xmesh, ymesh);
bc.left.a = 0; bc.left.b=1; bc.left.c=0; % left boundary condition
bc.right.a = 0; bc.right.b=1; bc.right.c=0; % right boundary condition
bc.bottom.a = 0; bc.bottom.b=1; bc.bottom.c=0; % bottom boundary condition
bc.top.a = 0; bc.top.b=1; bc.top.c=0; % top boundary condition

% Set the time vector and solve the PDE
tspan = linspace(0,10,50);
sol = pdepe(0,@pdex4pde,@pdex4ic,@pdex4bc,msh,tspan);

% Extract the solution vector from the solution matrix
u = sol(:,:,1);

% Plot the results
surf(x,y,u);
xlabel('x');
ylabel('y');
zlabel('u');
930 chars
31 lines

main.m
% Define the PDE
function [c,f,s] = pdex4pde(x,y,t,u,DuDx,DuDy)
c = 1;
f = D*[DuDx,DuDy];
s = 0;

% Define the initial conditions
function u0 = pdex4ic(x,y)
u0 = sin(pi*x).*sin(pi*y);

% Define the boundary conditions
function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)
pl = ul - 0;
ql = 0;
pr = ur - 0;
qr = 0;
308 chars
17 lines