Thank you Gianluca.
I am trying this variant:
K1 = i;
V = Tan[x];
P = r^2/(Cos[y])^5
eq = K1*r^3*D[u[r, y], {r, 3}] + r^2*D[u[r, y], {r, 2}] + r*K1*D[V^2*u[r, y], {r, 1}] + V*r*K1*D[u[r, y], {r, 1}] + V*r*K1*D[V*u[r,y], {r,1}]+ V^2*r*K1*D[u[r, y], {r, 1}] +
V^3*K1*D[u[r, y], {r, 1}] == 2*u[r,y]*P
sol = u[r, y] /.
NDSolve[{eq,
u[0, y] == 1,
Derivative[1, 0][u][0, y] == 0,
Derivative[2, 0][u][0, y] == 10,
Derivative[3, 0][u][0, y] == 0,
u[r, 0] == 1,
Derivative[0, 1][u][r, 0] == 0,
Derivative[0, 2][u][r, 0] == 10,
Derivative[0, 3][u][0, y] == 0},
u, {r, 0, 3}, {y, 0, 3},
MaxSteps -> Infinity, PrecisionGoal -> 1,
AccuracyGoal -> 1,
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 32, "MaxPoints" -> 32, "DifferenceOrder" -> 2},
Method -> {"Adams", "MaxDifferenceOrder" -> 1}}] //
Plot3D[sol, {r, 0, 3}, {y, 0, 3}, AxesLabel -> Automatic]
However, I get that
NDSolve::dsvar: 0.00021449999999999998` cannot be used as a variable.