solution 2D POISSON equation

Alhamdulillah...
I hope Allah SWT keep me from USELESS knowledge, keep me from VAIN activity, keep me from NO SOLEMN feeling, keep me from unsatisfied desire, and keep me from unanswered prayer. amin

then, I would like to give enormous thank to Mr Sungkono M.Si who forces me to study computation physics and always inspire me in using MATLAB.

now we try to get solution of 2D Poisson equation using finite differences and we will visualise the solution. the governing equation is given by

we use regular grid of 5 X 5 points and the model size 50 X 50 km. A finite-difference representation of the Poisson equation in 2D can be derived

then, we will get desain of node that we need to find the solution in a node. see this figure,

then the general solution

note: you have to change general solution value indices from geometrical indices(i and j) to global index of unknown k.

where Ny is vertical resolution.

L=zeros(25,25);
p=50/4;
q=50/4;
for k=1:5
    L(k,k)=1;
end
for k=7:19
    n=5;
    L(k,k)=(-2/p)+(-2/q);
    L(k,k-1)=1/q;
    L(k,k+1)=1/p;
    L(k,k-n)=1/p;
    L(k,k+n)=1/q;
end
for k=21:25
    L(k,k)=1;
end
for k=1:3
        n=5;
        L(k*n+1,:)=0;
        L(k*n,:)=0;
end
for k=1:25
    L(k,k)=(-2/p)+(-2/q);
end
R=zeros(25,1);
for k=7:19
    R(k,1)=1;
end
for k=1:4
    n=5;
    R(k*n+1)=0;
    R(k*n)=0;
end
L
R
x=L\R
for i=1:1:5
  for j=1:1:5
    k       =   (j-1)*5+i;
    phi(i,j) =   x(k);
  end
end
phi

% Making vectors for nodal points positions
x=1:5; % Horizontal
y=1:5; % Vertical

% Making new figure
figure(1);
% Plotting solution
surf(x,y,phi);
light;
shading interp;
colorbar;
lighting phong;
title('Solution of 2D Poisson equation')
xlabel('x, km')
ylabel('y, km')
zlabel('Gravity potential, J/kg')