%% heat equation approximate solution
% u_t = u_xx
% on x \in [0, 1) with periodic boundary conditions
% and initial condition of u(x, 0) = sin(2*pi*x)

N = 20;
h = 1/N;
x = 0:h:(1-h);   % ! why not 1?
x = x';

e = ones(N, 1);

% build a sparse matrix
L = spdiags([e -2*e e], [-1 0 1], N, N);
L(1, N) = 1;
L(N, 1) = 1;

% full(L)   % to see it

L = (1/h^2) * L;

% IC
u0 = sin(2*pi*x);

%% Static test of 2nd-deriv approximation
uxx_exact = -4*pi^2 * u0;
uxx = L*u0;
err = uxx_exact - uxx;
[uxx_exact  uxx  err]
lw = 'linewidth';
figure(1); clf;
plot(x, uxx_exact, 'ro-', lw, 3);
hold on
plot(x, uxx, 'bx-', lw, 2)
legend('exact', 'approx')
xlabel('x');
ylabel('u_{xx}')


%% time-step size
k = 0.25*h^2;   % does this value matter?

Tf = 0.1
numsteps = ceil(Tf / k)
k = Tf/numsteps

u = u0;

%% loop to solve in time
for n = 1:numsteps
  unew = u + k*(L*u);

  if (mod(n, 10) == 0 || n == numsteps)
    % draw every 10th step
    figure(2); clf;
    plot(x, unew, 'kx-', lw, 3)
    hold on
    plot(x, u0, 'r--', lw, 3)
    xlabel('x'); ylabel('u')
    title(['t = ' num2str(n*k)]);
    drawnow
  end

  u = unew;
end
