% Two discretisations that share a matrix.
%
%   1. the three-point scheme:      (u_{i-1} - 2u_i + u_{i+1}) / h^2 = -f_i
%   2. the compact fourth-order scheme, due to Numerov:
%
%        (u_{i-1} - 2u_i + u_{i+1}) / h^2 = -(f_{i-1} + 10 f_i + f_{i+1}) / 12
%
% The left-hand sides are identical. All that changes is that the source is
% averaged over the same three points with weights 1/12, 10/12, 1/12, which is
% exactly the combination that cancels the h^2 term of the truncation error.
% The matrix is still tridiag(-1, 2, -1), the solve is still O(N), and the
% method is fourth-order accurate: the extra two orders cost one weighted
% average of numbers we already had.
%
% The catch is that the scheme wants f at the nodes and nothing else, so it
% needs f to be smooth as well: the cancellation is a statement about the
% Taylor expansion of f, not just of u.

function names = method_names()
names = {'3-point differences', 'compact 4th order (Numerov)'};
end

function [x, u] = solve1d(mid, N, f, ua, ub)
h = 1 / N;
x = (0:N)' / N;
fx = f(x);

if mid == 1
    r = h^2 * fx(2:N);
else
    r = h^2 * (fx(1:N-1) + 10*fx(2:N) + fx(3:N+1)) / 12;
end
r(1) = r(1) + ua;
r(end) = r(end) + ub;

u = [ua; thomas(N-1, r); ub];
end

function w = thomas(n, d)
c = zeros(n, 1);
e = zeros(n, 1);
b = 2;
c(1) = -1 / b;
e(1) = d(1) / b;
for i = 2:n
    b = 2 + c(i-1);
    c(i) = -1 / b;
    e(i) = (d(i) + e(i-1)) / b;
end
w = zeros(n, 1);
w(n) = e(n);
for i = n-1:-1:1
    w(i) = e(i) - c(i) * w(i+1);
end
end