function temple3043_heateqn
%TEMPLE3043_HEATEQN
%   Second order accurate numerical Appromation of the 2D heat
%   equation, based on the Crank-Nicolson time stepping scheme.
%   The method requies a linear system to be solved in each
%   time step, which is here done by Backslash.
%   The purpose of this code is to serve as a wrapper
%   application for the investigation of linear solvers. To do
%   do, replace the Backslash solve by other solvers, such as:
%   (1) precomputed direct methods (LU, Cholesky, sparse elim.)
%   (2) relaxation schemes (Jacobi, Gauss-Seidel, SOR)
%   (3) Krylov methods (Conjugate gradients)
%   Furthermore, other techniques can be investigated, such as
%   (4) preconditioning
%   (5) multigrid
%   (6) smart iteration starting vectors
%   (7) subspace recycling.
%
%   (C) 2011/10/23 by Benjamin Seibold

% Initialization
tf = 1e-2; % final time
n = 500; % number of cells in each dimension
h = 1/n; % distance between cell centers
dt = 1e-1*h; % desired time step
nt = ceil(tf/dt); dt = tf/nt; % actual time step
x = h/2:h:1-h/2; % vector of cell center positions
[X,Y] = meshgrid(x); % mesh for plotting
U = exp(-60*(X-.2).^2-15*(Y-.6).^2)+...
    exp(-15*(X-.6).^2-40*(Y-.7).^2)+...
    exp(-200*(X-.8).^2-20*(Y-.4).^2)+...
    exp(-30*(X-.5).^2-200*(Y-.2).^2); % initial conditions
u = reshape(U',n^2,1); % reshape 2D solution into column vector

%System matrix setup
e = ones(n,1); % vector of ones
K1D = spdiags(e*[1 -2 1],-1:1,n,n)/h^2; % 1D discrete Laplacian matrix
I1D = speye(n); % 1D identity matrix
K2D = kron(K1D,I1D)+kron(I1D,K1D); %2D discrete Laplacian
I2D = speye(n^2); % 2D identity matrix
A = I2D-(dt/2)*K2D; % matrices used for Crank-Nicolson
B = I2D+(dt/2)*K2D; % time-stepping

% Time loop
cost = 0;
for i = 1:nt
    % computation
    tic
    b = B*u; % compute b = B*u my matrix multiplication
    u = A\b; % compute u = inv(A)*b by solving the linear system with A
    cost = cost+toc;
    % visualization
    U = reshape(u,n,n)'; % reshape solution into 2d array for plotting
    clf, surf(X,Y,U,'EdgeAlpha',.1) % plot solution
    axis([0 1 0 1 0 1]), caxis([0 1]) % axes for solution and color
    xlabel('x'), ylabel('y'), zlabel('u(x,y,t)') % axes labels
    title(sprintf('Heat equation at t=%0.3f',i*dt))
    view(-20,50), drawnow
end

 % Report run time of solution
fprintf('Time spent on computation: %0.3f seconds\n',cost)
