clear; clc; close all;

n = 5; V0 = 1;

x = linspace(0, 1, 5); y = linspace(0, 1, 5);
[X, Y] = meshgrid(x, y);

figure; hold on; grid on;
plot(X, Y, 'k.', 'markersize', 24);
for i = 1 : length(x)
    for j = 1 : length(y)
        text(X(n * (i - 1) + j), Y(n * (i - 1) + j) + 0.05, sprintf('V%d', n * (i - 1) + j));
    end
end

% boundary condition
phi = zeros(1, length(x) * length(y));
phi(5 : 5 : 25) = V0;

A = [4 -1 0 -1 0 0;
    -1 4 -1 0 -1 0;
    0 -1 4 0 0 -1;
    -2 0 0 4 -1 0;
    0 -2 0 -1 4 -1;
    0 0 -2 0 -1 4];
b = [0; 0; V0; 0; 0; V0];

% inverse of the matrix
v = A \ b

% generate the solution matrix
phi([7 8 9]) = v(1 : 3);
phi([17 18 19]) = phi([7 8 9]); % by symmetry
phi([12 13 14]) = v(4 : 6);

phi = reshape(phi, 5, 5);
for i = 1 : length(y)
    for j = 1 : length(x)
        h = text(X(n * (i - 1) + j), Y(n * (i - 1) + j) - 0.05, sprintf('%7.4f', phi(j, i)));
        set(h, 'color', 'r');
    end
end

contour(X, Y, phi, 10);