% Second-order finite differences for  -u'' = f  on (0, 1),  u(0) = ua, u(1) = ub.
%
% Replace u''(x_i) by the three-point second difference
%
%     u''(x_i) ~ ( u(i-1) - 2 u(i) + u(i+1) ) / h^2 ,
%
% which is exact for quadratics and wrong by (h^2/12) u''''(x_i) for anything
% else. Collecting the N-1 interior equations gives a tridiagonal system whose
% matrix is (1/h^2) * tridiag(-1, 2, -1); the two known boundary values move to
% the right-hand side.

function names = method_names()
names = {'3-point differences'};
end

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

r = h^2 * f(x(2:N));      % scale by h^2 so the matrix is tridiag(-1, 2, -1)
r(1) = r(1) + ua;         % the neighbour at x = 0 is known
r(end) = r(end) + ub;     % and so is the one at x = 1

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

function w = thomas(n, d)
% Solve tridiag(-1, 2, -1) w = d by Gaussian elimination without pivoting,
% which is safe here because the matrix is symmetric positive definite and
% diagonally dominant. Two sweeps, O(n) work and O(n) storage: the whole
% reason a one-dimensional problem is cheap at any resolution.
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);       % = 2 - (-1)*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