function poisson2d_convergence()
% POISSON2D_CONVERGENCE  Convergence test for the 2D Poisson equation
%   Solves  -Laplacian(u) = f  on [0,1]x[0,1]
%   with homogeneous Dirichlet BCs using 2nd-order finite differences.
%
%   Manufactured solution: u(x,y) = sin(2*pi*x) * sin(2*pi*y)
%   so that f(x,y) = 8*pi^2 * sin(2*pi*x) * sin(2*pi*y).
%
%   Runs on grids N = 16,32,64,128,256 and reports L2/Linf errors
%   and convergence rates. Expects 2nd-order (rate ~ 2).

    % --- Setup -----------------------------------------------------------
    Nvec = [16 32 64];
    errL2   = zeros(size(Nvec));
    errLinf = zeros(size(Nvec));

    % Exact solution and RHS
    uExact = @(x,y) sin(2*pi*x) .* sin(2*pi*y);
    fRHS   = @(x,y) 8*pi^2 * sin(2*pi*x) .* sin(2*pi*y);

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

    % --- Solve on each grid ----------------------------------------------
    for k = 1:numel(Nvec)
        N  = Nvec(k);
        h  = 1 / (N + 1);              % grid spacing
        x  = h * (1:N);                % interior nodes
        y  = h * (1:N);
        [X, Y] = meshgrid(x, y);       % N-by-N interior grid

        % Assemble the 2D Laplacian via Kronecker products
        e  = ones(N, 1);
        T  = spdiags([e -2*e e], -1:1, N, N) / h^2;   % 1D Laplacian
        I  = speye(N);
        A  = -(kron(I, T) + kron(T, I));   % -Laplacian (N^2 x N^2)

        % Right-hand side (column-major ordering matches meshgrid)
        F  = fRHS(X, Y);
        b  = F(:);

        % Solve the linear system
        u  = A \ b;
        U  = reshape(u, N, N);         % numerical solution on interior

        % Exact solution on interior
        Ue = uExact(X, Y);

        % Errors
        errLinf(k) = max(abs(U(:) - Ue(:)));
        errL2(k)   = h * norm(U(:) - Ue(:));   % scaled discrete L2

        % --- Animate: show each grid result ------------------------------
        subplot(1, 2, 1);
        % Embed interior in full grid (with zero BCs on boundary)
        xf = [0, x, 1];
        yf = [0, y, 1];
        [Xf, Yf] = meshgrid(xf, yf);
        Ufull = zeros(N+2, N+2);
        Ufull(2:end-1, 2:end-1) = U;
        surf(Xf, Yf, Ufull, 'EdgeColor', 'none');
        colorbar;
        title(sprintf('Numerical solution  (N = %d)', N));
        xlabel('x'); ylabel('y'); zlabel('u');
        view(135, 30);

        subplot(1, 2, 2);
        Efull = zeros(N+2, N+2);
        Efull(2:end-1, 2:end-1) = abs(U - Ue);
        surf(Xf, Yf, Efull, 'EdgeColor', 'none');
        colorbar;
        title(sprintf('|Error|  (N = %d,  max = %.2e)', N, errLinf(k)));
        xlabel('x'); ylabel('y'); zlabel('|error|');
        view(135, 30);

        drawnow;
        pause(1.0);
    end

    % --- Convergence rate plot -------------------------------------------
    hvec = 1 ./ (Nvec + 1);

    figure('Position', [100 100 600 500]);
    loglog(hvec, errLinf, 'bo-', 'LineWidth', 1.5, 'MarkerSize', 8); hold on;
    loglog(hvec, errL2,   'rs-', 'LineWidth', 1.5, 'MarkerSize', 8);
    loglog(hvec, hvec.^2 * (errLinf(1)/hvec(1)^2), 'k--', 'LineWidth', 1);
    legend('L_\infty error', 'L_2 error', 'O(h^2) ref', 'Location', 'southeast');
    xlabel('h'); ylabel('Error');
    title('Convergence of 2D Poisson Solver');
    set(gca, 'FontSize', 12);
    grid on;
    saveas(gcf, 'poisson2d_convergence.png');

    % --- Summary table ---------------------------------------------------
    rateLinf = [NaN, log2(errLinf(1:end-1) ./ errLinf(2:end))];
    rateL2   = [NaN, log2(errL2(1:end-1)   ./ errL2(2:end))];

    fprintf('\n');
    fprintf('===================================================================\n');
    fprintf('  2D Poisson Convergence Test — u(x,y)=sin(2πx)sin(2πy)\n');
    fprintf('===================================================================\n');
    fprintf('%6s %10s %12s %8s %12s %8s\n', ...
            'N', 'h', 'L_inf err', 'rate', 'L2 err', 'rate');
    fprintf('-------------------------------------------------------------------\n');
    for k = 1:numel(Nvec)
        if k == 1
            fprintf('%6d %10.4e %12.4e %8s %12.4e %8s\n', ...
                    Nvec(k), hvec(k), errLinf(k), '--', errL2(k), '--');
        else
            fprintf('%6d %10.4e %12.4e %8.4f %12.4e %8.4f\n', ...
                    Nvec(k), hvec(k), errLinf(k), rateLinf(k), errL2(k), rateL2(k));
        end
    end
    fprintf('===================================================================\n');
    fprintf('Expected rate: 2.00  (2nd-order central differences)\n\n');
end
