potential in 2D space by finite difference method

% Define parameters
Nx = 50; % Number of grid points in x-direction
Ny = 50; % Number of grid points in y-direction
Lx = 1; % Length of domain in x-direction
Ly = 1; % Length of domain in y-direction
dx = Lx / (Nx - 1); % Grid spacing in x-direction
dy = Ly / (Ny - 1); % Grid spacing in y-direction

% Initialize potential grid
phi = zeros(Ny, Nx); % Initialize potential grid with zeros

% Boundary conditions
phi(1, :) = 0; % Boundary condition at y = 0
phi(Ny, :) = 0; % Boundary condition at y = Ly
phi(:, 1) = 0; % Boundary condition at x = 0
phi(:, Nx) = 1; % Boundary condition at x = Lx

% Finite difference loop
tolerance = 1e-6; % Tolerance for convergence
maxIterations = 1000; % Maximum number of iterations
iteration = 0; % Initialize iteration counter
while true
    oldPhi = phi; % Store old potential values
    
    % Update potential using finite difference method
    for i = 2:Nx-1
        for j = 2:Ny-1
            phi(j, i) = (phi(j+1, i) + phi(j-1, i) + phi(j, i+1) + phi(j, i-1)) / 4;
        end
    end
    
    % Check for convergence
    residual = max(max(abs(phi - oldPhi))); % Compute residual
    if residual < tolerance || iteration >= maxIterations
        break; % Exit loop if converged or reached maximum iterations
    end
    
    iteration = iteration + 1; % Increment iteration counter
end

% Display results
fprintf('Number of iterations: %d\n', iteration);
fprintf('Residual: %e\n', residual);

% Plot potential
[X, Y] = meshgrid(0:dx:Lx, 0:dy:Ly);
surf(X, Y, phi);
xlabel('x');
ylabel('y');
zlabel('Potential');
title('Potential in 2D Space');