Laplace

Exact solution and an F

{a,b}={0.,3};
{c,d}={0.,12};
{nx,ny}={204,206};
{Δx,Δy}={((b-a)/nx),((d-c)/ny)};
u[x_,y_]= Cos[π/2 x] Sin[y π/3]
f[x_,y_]= D[ u[x,y],{x,2}] +D[ u[x,y],{y,2}] 
UEx=Table[ 
u[x,y],
{y,c,d,Δy},
{x,a,b,Δx}
];
FMat=Table[ 
f[x,y],
{y,c,d,Δy},
{x,a,b,Δx}
];
F=Flatten[FMat];
{ListPlot3D[UEx],Plot3D[ u[x,y],{x,a,b},{y,c,d}]}

D[u[x,y],x]/.x->a
u[b,y]

Matrix

kFun[i_,j_]:=i+(nx+1)*(j -1)
Table[
kFun[i,j],
{j,1,ny+1},
{i,1,nx+1}
];
(* Comment *)
IntNodes=Flatten[
Table[
k=kFun[i,j];
{
{k,k}->-(2/((Δx)^2))-(2/((Δy)^2)),
{k,k+1}->(1/((Δx)^2)),
{k,k-1}->(1/((Δx)^2)),
{k,kFun[i,j+1]}->(1/((Δy)^2)),
{k,kFun[i,j-1]}->(1/((Δy)^2))

},
{i,2,nx},
{j,2,ny}
]
];
NeumanmNodes=Flatten[
Table[
k=kFun[1,j];
F[[k]]=0;
{{k,k}->1,{k,k+1}->-1},
{j,1,ny+1}
]
];
DirichletNodes=Flatten[
Table[
k=kFun[nx+1,j];
F[[k]]=0;
{k,k}->1,
{j,1,ny+1}
]
];
TopNodes=Flatten[
Table[
k=kFun[i,ny+1];
{
{k,k}->-(2/((Δx)^2))-(2/((Δy)^2)),
{k,k+1}->(1/((Δx)^2)),
{k,k-1}->(1/((Δx)^2)),
{k,kFun[i,ny+1]}->(1/((Δy)^2)),
{k,kFun[i,1]}->(1/((Δy)^2))

},
{i,2,nx}
]
];
BottomNodes=Flatten[
Table[
k=kFun[i,1];
{
{k,k}->-(2/((Δx)^2))-(2/((Δy)^2)),
{k,k+1}->(1/((Δx)^2)),
{k,k-1}->(1/((Δx)^2)),
{k,kFun[i,2]}->(1/((Δy)^2)),
{k,kFun[i,ny+1]}->(1/((Δy)^2))

},
{i,2,nx}
]
];
M=SparseArray[
Join[IntNodes,DirichletNodes, NeumanmNodes, TopNodes, BottomNodes]
]

Q=LinearSolve[M, Flatten[F]];
QMat=Partition[Q,nx+1];
{ListPlot3D[QMat],ListPlot3D[UEx]}
ListPlot3D[QMat-UEx]
ListPlot[ Flatten[ QMat-UEx]]