Your code should somehow look like this. Fill in the details.
nx = 11;
ny = 11;
Lx = 0.1;
Ly = 0.1;
dx = Lx/(nx-1);
dy = Ly/(ny-1);
for ix = 1:nx
for iy = 1:ny
T(ix,iy) = 600;
end
end
for ix = 1:nx
for iy = 1:ny
Y0((ix-1)*ny + iy ) = T(ix,iy);
end
end
[time,Y] = ode15s(@(time,Y)fun(time,Y,nx,ny,dx,dy),tspan,Y0)
function dYdt = fun(time,Y,nx,ny,dx,dy)
% Write solution vector Y in matrix
T = zeros(nx,ny);
for ix = 1:nx
for iy = 1:ny
T(ix,iy) = Y((ix-1)*ny + iy);
end
end
% Allocate workspace for the time derivatives in the grid points
dTdt = zeros(nx,ny);
% Set the dTdt expressions of your attached paper (Don't use function
% handles, but just the variables here !)
...
% Write back the dTdt matrix into a column vector to return it to ODE15S
for ix = 1:nx
for iy = 1:ny
dYdt((ix-1)*ny + iy) = dTdt(ix,iy);
end
end
end
