%% Poisson Equation Demo - Finite Difference Method
% Solves -Laplacian(u) = f on [0,1]x[0,1]
% with homogeneous Dirichlet boundary conditions
% Source term: f = 2*pi^2 * sin(pi*x) * sin(pi*y)
% Exact solution: u = sin(pi*x) * sin(pi*y)

clear; close all;

% Grid parameters
N = 50;                      % Number of interior points per dimension
h = 1/(N+1);                 % Grid spacing
x = linspace(0, 1, N+2);     % x coordinates including boundary
y = linspace(0, 1, N+2);     % y coordinates including boundary
[X, Y] = meshgrid(x, y);

% Source term
f = @(x,y) 2*pi^2 * sin(pi*x) .* sin(pi*y);

% Exact solution for comparison
u_exact = @(x,y) sin(pi*x) .* sin(pi*y);

% Build the coefficient matrix using 5-point stencil
% -u_{i-1,j} - u_{i+1,j} - u_{i,j-1} - u_{i,j+1} + 4*u_{i,j} = h^2 * f_{i,j}
e = ones(N, 1);
T = spdiags([-e, 4*e, -e], [-1, 0, 1], N, N);  % 1D Laplacian block
I = speye(N);
A = kron(I, T) + kron(spdiags([-e, -e], [-1, 1], N, N), I);

% Build the right-hand side vector
F = zeros(N, N);
for i = 1:N
    for j = 1:N
        F(i,j) = h^2 * f(x(i+1), y(j+1));
    end
end
b = F(:);

% Solve the linear system
u_interior = A \ b;

% Reshape solution to 2D grid
U = zeros(N+2, N+2);
U(2:N+1, 2:N+1) = reshape(u_interior, N, N);

% Compute exact solution on grid
U_exact = u_exact(X, Y);

% Compute error
error_max = max(max(abs(U - U_exact)));
error_L2 = h * norm(U(:) - U_exact(:));

%% Create figure
fig = figure('Position', [100, 100, 1200, 500]);

% Plot numerical solution
subplot(1,3,1);
surf(X, Y, U);
xlabel('x'); ylabel('y'); zlabel('u');
title('Numerical Solution');
shading interp;
colorbar;
view(45, 30);

% Plot exact solution
subplot(1,3,2);
surf(X, Y, U_exact);
xlabel('x'); ylabel('y'); zlabel('u');
title('Exact Solution');
shading interp;
colorbar;
view(45, 30);

% Plot error
subplot(1,3,3);
surf(X, Y, abs(U - U_exact));
xlabel('x'); ylabel('y'); zlabel('|error|');
title(sprintf('Error (max = %.2e)', error_max));
shading interp;
colorbar;
view(45, 30);

sgtitle(sprintf('Poisson Equation: -\\Delta u = f, N = %d, h = %.4f', N, h));

% Save figure
print(fig, 'poisson_result', '-dpng', '-r150');
fprintf('Saved poisson_result.png\n');
fprintf('Max error: %.6e\n', error_max);
fprintf('L2 error:  %.6e\n', error_L2);
