(* 
(((∂^2)u)/(∂(x^2)))+(((∂^2)u)/(∂(y^2)))=f(x,y);

u(0,y)=ϕ[y], u(1,y)=ψ(y);
u(x,0)=α[x], u(x,1)=β(x);

(Subscript[u, rx,ry])=u(rx Δx,ry Δy)


 *)

Plot3D[f[x,y],{x,0,1},{y,0,1},PlotRange->All]

ϕ[y_]:=0;
ψ[y_]:=0;
α[x_]:=0;
β[x_]:=0;
f[x_,y_]:=-Exp[-((((x-0.5))^2)+(((y-0.5))^2))/(2 ((0.1^2)))];
nx=30;
ny=30;
Δx=1.0/nx;
Δy=1.0/ny;
reg1=Table[Subscript[u,0,ry]->ϕ[ry Δy],{ry,1,ny-1}];
reg2=Table[Subscript[u,nx,ry]->ψ[ry Δy],{ry,1,ny-1}];
reg3=Table[Subscript[u,rx,0]->α[rx Δx],{rx,1,nx-1}];
reg4=Table[Subscript[u,rx,ny]->β[rx Δx],{rx,1,nx-1}];
reg=Join[reg1,reg2,reg3,reg4];

eqs=Flatten[Table[
f[rx Δx,ry Δy]==(Subscript[u,rx+1,ry]-2Subscript[u,rx,ry]+Subscript[u,rx-1,ry])/((Δx)^2)+(Subscript[u,rx,ry+1]-2Subscript[u,rx,ry]+Subscript[u,rx,ry-1])/((Δy)^2),{rx,1,nx-1},{ry,1,ny-1}]]/.reg;

vars=Flatten[Table[Subscript[u,rx,ry],{rx,1,nx-1},{ry,1,ny-1}]];

sol=NSolve[eqs,vars];

data={};
Do[
vx=rx Δx;
vy=ry Δy;
vu=Subscript[u,rx,ry]/.sol[[1]];
point={vx,vy,vu};
AppendTo[data,point];
,{rx,1,nx-1},{ry,1,ny-1}];
Do[
vx=0.0;
vy=ry Δy;
vu=ϕ[vy];
point={vx,vy,vu};
AppendTo[data,point];
vx=1.0;
vy=ry Δy;
vu=ψ[vy];
point={vx,vy,vu};
AppendTo[data,point];
,{ry,0,ny}];
Do[
vx=rx Δx;
vy=0.0;
vu=α[vx];
point={vx,vy,vu};
AppendTo[data,point];
vx=rx Δx;
vy=1.0;
vu=β[vx];
point={vx,vy,vu};
AppendTo[data,point];
,{rx,0,nx}];

ListPlot3D[data,PlotRange->All]

data