%Example Code by Philip Caplan
%https://de.mathworks.com/matlabcentral/profile/authors/5560510-philip-caplan

function solve_heat
clear;clc;close all;

% compute until which time?
tmax  = .1;

% make the grid
nx = 11;
ny = 11;
[x,y] = meshgrid(linspace(0,1,nx),linspace(0,1,ny));

% initial condition
u0 = sin(pi*x).*sin(pi*y);
u0 = reshape(u0,[nx*ny,1]);

du_dt = @(t,u) heat_rhs(x,y,t,u);

u = ode45(du_dt,[0,tmax],u0);

end

function du_dt = heat_rhs(x,y,t,u)
% --- Computes du/dt = kVu at interior points
% 
% Inputs:
%   x,y  : grid coordinates
%   t    : current time
%   u    : solution at previous time step
%
% Outputs:
%   dudt : du/dt at current time step

% get some sizes
nx = size(x,1);
ny = size(x,2);

% reshape since ode45 handles a vector
u = reshape(u,[nx,ny]);

% plot solution at current time step
surf(x,y,u);
caxis([-1,1]);
zlim([-1,1]);
colorbar;
drawnow;

% assume uniform grid
dx = x(1,2) -x(1,1);
dy = y(2,1) -y(1,1);

% [TREAT BOUNDARY CONDITIONS HERE]
% (not needed in this case for homogeneous Dirichlet bc's)

d2u_dx2 = zeros(nx,ny);
d2u_dy2 = zeros(nx,ny);

% second derivative operator
d2u_dx2(2:end-1,:) = (u(3:end,:) -2*u(2:end-1,:) +u(1:end-2,:))/dx^2;
d2u_dy2(:,2:end-1) = (u(:,3:end) -2*u(:,2:end-1) +u(:,1:end-2))/dy^2;

du_dt = d2u_dx2 +d2u_dy2;

% reshape back to vector for ode45
du_dt = reshape(du_dt,[nx*ny,1]);

end