% The leapfrog scheme for the wave equation u_tt = c^2 u_xx.
%
% Centred second differences in both time and space:
%
%   u(j, n+1) = 2 u(j, n) - u(j, n-1) + lam^2 (u(j+1, n) - 2 u(j, n) + u(j-1, n)),
%
% with lam = c dt / h the Courant number.  It is second-order accurate, it is
% stable exactly when lam <= 1, and it has a property that no other scheme in
% this tutorial shares: at lam = 1 it is not an approximation at all.
%
% u is the displacement at step 0 and uold at step -1, one entry per node.
% bc is 'fixed', in which case u has N+1 entries and the first and last are
% held at zero, or 'periodic', in which case it has N entries and node N is
% followed by node 1.  The scheme advances M steps and returns
%
%   y  M-by-1         the displacement at node jr after each step
%   U  numel(u)-by-K  the whole field after each of the first K steps
%
% The inner loop is written out in full, one node at a time; numbl compiles it,
% and in this form it is about six times faster than the vectorised slice.

function [y, U] = leapfrog(u, uold, lam, M, bc, jr, K)
u = u(:);
uold = uold(:);
n1 = numel(u);
l2 = lam ^ 2;
y = zeros(M, 1);
U = zeros(n1, K);
unew = zeros(n1, 1);
periodic = strcmp(bc, 'periodic');
for n = 1:M
    for j = 2:n1 - 1
        unew(j) = 2 * u(j) - uold(j) + l2 * (u(j + 1) - 2 * u(j) + u(j - 1));
    end
    if periodic
        unew(1) = 2 * u(1) - uold(1) + l2 * (u(2) - 2 * u(1) + u(n1));
        unew(n1) = 2 * u(n1) - uold(n1) + l2 * (u(1) - 2 * u(n1) + u(n1 - 1));
    else
        unew(1) = 0;
        unew(n1) = 0;
    end
    uold = u;
    u = unew;
    y(n) = u(jr);
    if n <= K
        U(:, n) = u;
    end
end
end