% The same five-point system, solved two ways.
%
% The matrix is A = kron(I,T) + kron(T,I) with T = tridiag(-1,2,-1), and the
% eigenvectors of T are known exactly: the sine vectors
%
%     v_k(j) = sin(pi j k / N),     T v_k = lambda_k v_k,
%     lambda_k = 2 (1 - cos(pi k / N)),   k = 1, ..., N-1.
%
% So the eigenvectors of A are the products v_k(i) v_l(j) and its eigenvalues
% are the sums lambda_k + lambda_l. Expanding the right-hand side in that basis,
% dividing by the eigenvalues and transforming back solves the system exactly,
% and the expansion is a discrete sine transform, which one FFT computes in
% O(N log N) per line. Total: O(N^2 log N) for N^2 unknowns, no matrix formed
% and nothing stored beyond the grid itself.
%
% This is not an iterative method and not an approximation. It returns the same
% vector the dense solve returns, to within roundoff, and it does so for grids
% the dense solve cannot represent.

function names = method_names()
names = {'5-point, dense solve', '5-point, sine transform'};
end

function [x, y, U] = solve2d(mid, N, g)
h = 1 / N;
x = (0:N)' / N;
y = x;
n = N - 1;

U = zeros(N+1, N+1);
U(1, :) = g(x', y(1));
U(N+1, :) = g(x', y(N+1));
U(:, 1) = g(x(1), y);
U(:, N+1) = g(x(N+1), y);

B = zeros(n, n);
B(1, :) = B(1, :) + U(1, 2:N);
B(n, :) = B(n, :) + U(N+1, 2:N);
B(:, 1) = B(:, 1) + U(2:N, 1);
B(:, n) = B(:, n) + U(2:N, N+1);

if mid == 1
    T = 2*eye(n) - diag(ones(n-1, 1), 1) - diag(ones(n-1, 1), -1);
    A = kron(eye(n), T) + kron(T, eye(n));
    U(2:N, 2:N) = reshape(A \ B(:), n, n);
else
    lam = 2 * (1 - cos(pi * (1:n)' / N));           % eigenvalues of T
    D = lam * ones(1, n) + ones(n, 1) * lam';       % eigenvalues of A
    Bh = lp_dst1(lp_dst1(B).').';                   % to the sine basis, both ways
    U(2:N, 2:N) = lp_dst1(lp_dst1(Bh ./ D).').';    % divide, and back again
end
end