# How to improve on the evaluation time

I want to solve a complicated set of 9 coupled PDEs in 2 dimensions on a Chebyshev lattice. As a simplified problem, I've started by writing the following code to solve one simple equation with uniform Dirichlet boundary conditions, but I need to make it as fast as possible. Any ideas on how to improve on the evaluation time?

Thanks a lot! Chris

(*getting the lattice*)

Ny = 30; Nx = 30;

CGLGrid[x0_, L_, n_Integer /; n > 1] :=
x0 + L/2 (1 - Cos[\[Pi] Range[0, n - 1]/(n - 1)]);
gridy = CGLGrid[-1, 2, Ny];
gridx = CGLGrid[-1, 2, Nx];

grid = Transpose[Partition[Tuples[{gridy, gridx}], Nx]];

(*getting the differentiation matrices*)

d[1, 0] =
NDSolveFiniteDifferenceDerivative[{1, 0}, {gridx, gridy},
"DifferenceOrder" -> {"Pseudospectral", "Pseudospectral"},
PeriodicInterpolation -> {False, False}];

d[0, 1] =
NDSolveFiniteDifferenceDerivative[{0, 1}, {gridx, gridy},
"DifferenceOrder" -> {"Pseudospectral", "Pseudospectral"},
PeriodicInterpolation -> {False, False}];

d[2, 0] =
NDSolveFiniteDifferenceDerivative[{2, 0}, {gridx, gridy},
"DifferenceOrder" -> {"Pseudospectral", "Pseudospectral"},
PeriodicInterpolation -> {False, False}];

d[0, 2] =
NDSolveFiniteDifferenceDerivative[{0, 2}, {gridx, gridy},
"DifferenceOrder" -> {"Pseudospectral", "Pseudospectral"},
PeriodicInterpolation -> {False, False}];

d[1, 1] =
NDSolveFiniteDifferenceDerivative[{1, 1}, {gridx, gridy},
"DifferenceOrder" -> {"Pseudospectral", "Pseudospectral"},
PeriodicInterpolation -> {False, False}];

(*discretising*)

Off[CompiledFunction::cfsa];

Eqnlist = {D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}] - 10 f[x, y]};
ff[x_, y_] := Sin[8 x (y - 1)];

Eqnlistloc =
Compile[{{Uval, _Real, {Nx, Ny}}},
Eqnlist /. {D[u[x, y], {y, 2}] -> d[0, 2][Uval],
D[u[x, y], {x, 2}] -> d[2, 0][Uval],
f[x, y] -> Partition[Apply[ff, Flatten[grid, 1], 1], Nx]}];

Eqnlistdis[Uval_] :=
Join[Flatten[
Transpose[
Take[Transpose[Take[Eqnlistloc[Uval][[1]], {2, Ny - 1}]], {2,
Nx - 1}]]]];

(*Newton method*)

Uval = ConstantArray[0., {Nx, Ny}];
For[i = 1, i < 40, i++;
M1 = Eqnlistdis[Uval];
dEqnlistdis = {};
For[ n = 1, n < Ny - 1, n++;
For[ m = 1, m < Nx - 1, m++;
Uval[[n]][[m]] = Uval[[n]][[m]] + 1/1000;
M2 = Eqnlistdis[Uval];
dEqnlistdis = Join[dEqnlistdis, {1000 (M2 - M1)}];
Uval[[n]][[m]] = Uval[[n]][[m]] - 1/1000;]];
dUval = Transpose[
Join[{ConstantArray[0., Ny]},
Transpose[
Join[{ConstantArray[0., Nx - 2]},
Partition[LinearSolve[dEqnlistdis, -M1],
Nx - 2], {ConstantArray[0., Nx - 2]}]], {ConstantArray[0.,
Ny]}]];
Uval = Uval + dUval;]
`
-