% Two methods on the same tridiagonal solver:
%
%   1. the three-point scheme of the previous chapter;
%   2. the same scheme run at N and at 2N, combined by Richardson
%      extrapolation to cancel the leading h^2 term of the error.
%
% If the error has an asymptotic expansion  u_h = u + c h^2 + O(h^4),  then
%
%     (4 u_{h/2} - u_h) / 3 = u + O(h^4)
%
% at the points the two meshes share. Nothing about the discretisation
% changes; the fourth-order accuracy comes out of two second-order solves and
% one subtraction. The expansion is what makes it work, and the expansion
% needs u to be smooth: on the low-regularity solution there is no h^2 term to
% cancel and the extrapolation gains nothing.

function names = method_names()
names = {'3-point differences', 'Richardson from N and 2N'};
end

function [x, u] = solve1d(mid, N, f, ua, ub)
if mid == 1
    [x, u] = fd2(N, f, ua, ub);
else
    [x, uc] = fd2(N, f, ua, ub);
    [~, uf] = fd2(2*N, f, ua, ub);
    u = (4 * uf(1:2:end) - uc) / 3;
end
end

function [x, u] = fd2(N, f, ua, ub)
h = 1 / N;
x = (0:N)' / N;
r = h^2 * f(x(2:N));
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