% The five-point Laplacian on the unit square, assembled as a matrix and
% solved by Gaussian elimination.
%
% On a uniform grid of N intervals in each direction, replace
%
%     u_xx + u_yy  by  ( u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1) - 4 u(i,j) ) / h^2
%
% at the (N-1)^2 interior points. Scaling by h^2 leaves the matrix
%
%     A = kron(I, T) + kron(T, I),      T = tridiag(-1, 2, -1),
%
% which is what the vec of the interior array is multiplied by. Boundary
% neighbours are known, so they move to the right-hand side.
%
% This is the honest, naive version: the matrix is formed in full and handed to
% backslash. It is also the version that stops working first. The unknown count
% is (N-1)^2, so the dense matrix has (N-1)^4 entries and the elimination costs
% (N-1)^6 / 3 operations: N = 32 needs 7 MB and a fraction of a second, N = 64
% needs 450 MB, and N = 128 is out of reach on any machine.

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

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

% boundary values, and zeros in the interior
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);

% right-hand side: h^2 f (zero here, since we are solving Laplace's equation)
% plus the known neighbours of the points next to the boundary
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);

T = 2*eye(n) - diag(ones(n-1, 1), 1) - diag(ones(n-1, 1), -1);
I = eye(n);
A = kron(I, T) + kron(T, I);

W = A \ B(:);
U(2:N, 2:N) = reshape(W, n, n);
end