How to write 2D finite volume solver in a functional approach?

Question

I'm solving 2D poisson problem of a finite volume problem using TDMA.

The solution of the problem is given:

\[Phi][x_, y_] := x^2 y + x y^2 + 1

The definitions are the next:

(*Input parameters*)
nx = 12; x0 = -1.0; xl = 1.0;
ny = 10; y0 = -1.0; yl = 1.0;
\[CapitalGamma] = 1.;

(*Meshing parameters*)
dx = (xl - x0)/nx;
dy = (yl - y0)/ny;
dz = 1.0;
dv = dx dy dz;
x = Range[x0, xl, dx];
y = Range[y0, yl, dy];
xc = (x + RotateRight[x])[[2 ;;]]/2.;
yc = (y + RotateRight[y])[[2 ;;]]/2.;

(*Source terms*)
g0y = yc - yc^2 + 1;
gly = yc + yc^2 + 1;
g0x = xc - xc^2 + 1;
glx = xc + xc^2 + 1;

(*Surface elements*)
Se = Outer[dy*dz &, x[[2 ;;]], yc];
Sw = Outer[dy*dz &, x[[;; -2]], yc];
Sn = Outer[dx*dz &, xc, y[[2 ;;]]];
Ss = Outer[dx*dz &, xc, y[[;; -2]]];

(*Coeffincients *)
aE = \[CapitalGamma]*Se/dx;
aW = \[CapitalGamma]*Sw/dx;
aN = \[CapitalGamma]*Sn/dy;
aS = \[CapitalGamma]*Ss/dy;
aP = aE + aW + aN + aS;
sP = Outer[(2 #1 + 2 #2)*dv &, xc, yc];

(*Boundary values*)
aP[[All, -1]] = aP[[All, -1]] + aN[[All, -1]];
sP[[All, -1]] = sP[[All, -1]] + 2 aN[[All, -1]]*glx;
aN[[All, -1]] = 0;

aP[[All, 1]] = aP[[All, 1]] + aS[[All, 1]];
sP[[All, 1]] = sP[[All, 1]] + 2 aS[[All, 1]]*g0x;
aS[[All, 1]] = 0;



aP[[-1, All]] = aP[[-1, All]] + aE[[-1, All]];
sP[[-1, All]] = sP[[-1, All]] + 2 aE[[-1, All]]*gly;
aE[[-1, All]] = 0;

aP[[1, All]] = aP[[1, All]] + aW[[1, All]];
sP[[1, All]] = sP[[1, All]] + 2 aW[[1, All]]*g0y;
aW[[1, All]] = 0;

(*Gauss TDMA function*)
tDMA[a_, b_, c_, d_] := Module[{m},
  m = SparseArray[{Band[{1, 1}] -> b, Band[{2, 1}] -> a[[2 ;;]], 
     Band[{1, 2}] -> c[[;; -2]]}, {Length[b], Length[b]}];
  LinearSolve[m, d, "Method" -> "Krylov"]
  ]

As using 2D i have to map my solution. Finite volumes book give me the algoritms by using loops as follows:

(*Gauss TDMA by lines for the x line*)
For[k = 0, k < 100, k++,
 sP2 = sP + aN*RotateLeft[Tb, {0, 1}] + aS*RotateRight[Tb, {0, 1}];
 For[j = 1, j <= ny, j++,
  Tb[[All, j]] = 
   tDMA[-aW[[All, j]], aP[[All, j]], -aE[[All, j]], sP2[[All, j]]]
  ]
 ]

(*Gauss TDMA by lines for the y line*)
For[k = 0, k < 100, k++,
 sP2 = aE*RotateLeft[Tb, {1, 0}] + aW*RotateRight[Tb, {1, 0}] + sP;
 For[i = 1, i <= nx, i++,
  Tb[[i, All]] = 
   tDMA[-aS[[i, All]], aP[[i, All]], -aN[[i, All]], sP2[[i, All]]]
  ]
 ]

I think this is not efficient if I use a large mesh. So I would like to write the loops in a functional approach.

I tried somthin like this

Tb = Outer[1 &, xc, yc];
Tb = Nest[
   MapThread[
     tDMA, {-aS, aP, -aN, 
      aE*RotateLeft[#, {1, 0}] + aW*RotateRight[#, {1, 0}] + sP}] &, 
   Tb, 20];

But it only works for the y line, not for the x line

For the xlines i tried trasposing all the aP,aE,aW,aN,aS, and sP arrays but also I think is not the most efficient solution.

(*Line X*)

s2 = Transpose[
    Nest[MapThread[
       tDMA, {-aWT, aPT, -aET, 
        aNT*RotateLeft[#, {1, 0}] + aST*RotateRight[#, {1, 0}] + 
         sPT}] &, Transpose[Tb], 20]]; // AbsoluteTiming

Where aWT, aPT, aET, aNT, aST, sPT are transposed aP,aE,aW,aN,aS, and sP.