# Solve PDE by using Spatial Discretization for two dimensional independent spatial coordinates

I have a code given by

beta = 4;
l = 0;
mu = 0;
q = 1/2;
M = 3/4;

V = (r^2 - 2 M r + q)/(r^6 (r^4 +beta q)^2) (r^8 (r^2 l (l + 1) + mu^2 r^4 +
2 M r - 2 q) +beta q r^4 (6 r^2 + mu^2 r^4 - 12 M r + 6 q) +
beta^2 q^2 (2 r^2 - r^2 l (l + 1) - 6 M r + 4 q));
X = r + ((2 M^2 - q) ArcTan[(-M + r)/Sqrt[Abs[-M^2 + q]]])/Sqrt[
Abs[-M^2 + q]] + M Log[q - 2 M r + r^2];

fV[z0_?NumericQ] :=With[{z = SetPrecision[z0, 100 + 1]},
If[Abs[z] <= 12,Re[V /. FindRoot[X == z, {r, 10000001/10000000},
MaxIterations -> 10000, WorkingPrecision -> 100]], 0]]

xgrid =  Join[ Range[-1, 44/3] 3/2];
ygrid =  Join[Range[0, 44/3] 3/2];
sol1 = First[NDSolve[{-4*D[S[x, y], x, y] == fV[1/2 (y - x)]*S[x, y],
S[x, 3/2] == 1, S[-3/2, y] == Exp[-.25 (y - 3/2)^2]},
S, {x, -3/2, 21}, {y, 0, 21},Method -> {"MethodOfLines",
"SpatialDiscretization" -> { "TensorProductGrid",
"Coordinates" -> {xgrid, ygrid}}}, AccuracyGoal -> 5]]


If you run the code, you can find an error at the end which Mathematica does not accept xgrid and ygrid as the coordinate lists for each spatial dimension x and y. How can I tell to Mathematica to take xgrid and ygrid as the coordinate lists? Actually, a similar example is presented in Mathematica help center with an additional temporal coordinate for Burgers's equation. I have follwed it in my code, but I do not know why Mathematica give error in this case. Thanks a lot for your help.

• The method of lines can be applied to computations involving one temporal variable and one or more spatial variables. "Coordinates" can be used only for the spatial variable. So, you can specify xgrid or ygrid, but not both. Commented Aug 24, 2018 at 23:00

sol1 = First[NDSolve[{-4*D[S[x, y], x, y] == fV[1/2 (y - x)]*S[x, y],
`