function temple8023_heateqn2d
%TEMPLE8023_HEATEQN2D
%    Solves the 2d heat equation u_t=u_xx+u_yy with homogeneous
%    Dirichlet b.c. by finite differences.
%    Time-stepping is done via the Crank-Nicolson method,
%    and the time step is chosen proportional to the grid spacing.
%    The true solution is provided, and the error at the final
%    time is computed, as is the CPU time to compute.
%    By doubling nt, the scheme's order and cost scaling can be
%    tested.

% 10/2019 by Benjamin Seibold
% Feel free to modify for teaching and learning.

% Problem parameters
tf = 0.01; % final time
nt = 4; % number of time steps
cfl = .5; % propertionality between time step and grid spacing (k/h)

% Iniatialization
k = tf/nt; % time step
h = k/cfl; % desired mesh resolution
nx = ceil(1/h-1); % number of interior grid points
h = 1/(nx+1); % grid spacing
x = (h:h:1-h)'; % 1d spatial grid
[X,Y] = meshgrid(x); % 2d spatial grid
u = @(x,y,t) sin(4*pi*x).*sin(2*pi*y)*exp(-((4*pi)^2+(2*pi)^2)*t)+...
    sin(6*pi*x).*sin(4*pi*y)*exp(-((6*pi)^2+(4*pi)^2)*t);
U = reshape(u(X,Y,0),[],1); % initial conditions

% Matrix generation
K1D = spdiags(ones(nx,1)*[1 -2 1],-1:1,nx,nx)/h^2; % 1d Poisson matrix
I1D = speye(size(K1D)); % 1d identity matrix
K2D = kron(K1D,I1D)+kron(I1D,K1D); % 2d Poisson matrix
I2D = speye(size(K2D)); % 2d identity matrix
A = I2D-k/2*K2D; % LHS matrix for Crank-Nicolson
B = I2D+k/2*K2D; % RHS matrix for Crank-Nicolson

% Time loop
cpu = 0; % initialize CPU time
for j = 0:nt
    t = j*k; % current time
    if j>0 % update step
        tic
        U = A\(B*U);
        cpu = cpu+toc;
    end
    u_true = u(X,Y,t); % evaluate true solution
    u_max = max(max(abs(u_true))); % solution value upper bound
    ax = [0,1,0,1,[-1,1]*u_max]; % plot axes
    clf
    subplot(1,2,1)
    mesh(x,x,u_true), axis(ax)
    title(sprintf('True solution at t=%0.3f',t))
    subplot(1,2,2)
    mesh(x,x,reshape(U,nx,nx)), axis(ax)
    title(sprintf('Approximation at t=%0.3f',t))
    pause(1e-2)
end
er = max(max(abs(reshape(U,nx,nx)-u(X,Y,t)))); % error in max norm
fprintf('error = %0.2e\n',er)
fprintf('compute time = %0.1fs\n',cpu)
