solution of 2D poisson equation(no.2, chapter3, Taras Gerya)

thanks to Mr.Sungkono, Boss Asdi Presetyo, Boss Jordan Eko Setiyawan who have shared at whole night.
clear all;
clc;
L=zeros(1271,1271);
px=1000/30;
py=1500/40;
for i=43:1229
    n=41;
    L(i,i-1)=1/py;
    L(i,i+1)=1/px;
    L(i,i-n)=1/px;
    L(i,i+n)=1/py;
end
for i=1:29
        n=41;
        L(i*n+1,:)=0;
        L(i*n,:)=0;
end
for i=1:1271
    L(i,i)=(-2/px)+(-2/py);
end
R=zeros(1271,1);
for i=43:1228
    R(i,1)=1;
end
for i=1:29
    n=41;
    R(i*n+1)=0;
    R(i*n)=0;
end
L
R
x=L\R
for i=1:1:41
  for j=1:1:31
    k       =   (j-1)*41+i;
    phi(i,j) =   x(k);
  end
end
phi

% Making vectors for nodal points positions
x=1:31; % Horizontal
y=1:41; % 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')

I belive that what I write above have any lack. I sincerely receive any comment to make it better. thanks