How to specify the initial conditions with ParametricNDSolve

Question

I am still working on this problem. I am trying to use shooting method since I do not know the correct value for the initial derivative. I am following this answer. Now the code looks like this (I have removed the coupling terms to see if I could arrive to a solution)

ClearAll["Global`*"]
m = 1;
r[u_, v_] := 2 m + 2 m ProductLog[-((u v)/E)];
F[u_, v_] := (32 m^3)/r[u, v] Exp[-(r[u, v]/(2 m))];
P1[u_, v_, l_] = (D[F[u, v], u] D[F[u, v], v])/F[u, v]^2 + (
   D[r[u, v], v] D[F[u, v], u])/(r[u, v] F[u, v]) - (
   D[r[u, v], u] D[r[u, v], v])/r[u, v]^2 - D[F[u, v], u, v]/
   F[u, v] - 3 D[r[u, v], u, v]/r[u, v] - 
   F[u, v]/(4 r[u, v]^2) l (l + 1);
Q1[u_, v_, l_] = -(F[u, v]/(4 r[u, v]^2)) l (l + 1) + (
   D[F[u, v], u] D[F[u, v], v])/F[u, v]^2 + (
   D[F[u, v], v] D[r[u, v], u])/(F[u, v] r[u, v]) - (
   D[r[u, v], u] D[r[u, v], v])/r[u, v]^2 - D[F[u, v], u, v]/
   F[u, v] - 3 D[r[u, v], u, v]/r[u, v];
P2[u_, v_] = -((D[F[u, v], u] D[r[u, v], v])/(F[u, v] r[u, v])) + 
   D[r[u, v], u]^2/r[u, v]^2 + D[r[u, v], {u, 2}]/r[u, v];
Q2[u_, v_] = -((D[F[u, v], v] D[r[u, v], u])/(F[u, v] r[u, v])) + 
   D[r[u, v], v]^2/r[u, v]^2 + D[r[u, v], {v, 2}]/r[u, v];
V[u_, v_] = (Q2[u, v] - P2[u, v])/2 + 
   1/2 ((D[F[u, v], u] D[r[u, v], v] - D[F[u, v], v] D[r[u, v], u])/(
     F[u, v] r[u, v]));
W[u_, v_, 
   l_] = ((P1[u, v, l] + Q1[u, v, l])/2 + (Q2[u, v] + P2[u, v])/
     2) + (-3/4 (D[F[u, v], v] D[F[u, v], u])/F[u, v]^2 - 
     1/2 (D[r[u, v], u] D[F[u, v], v])/(r[u, v] F[u, v]) - 
     1/2 (D[r[u, v], v] D[F[u, v], u])/(r[u, v] F[u, v]) + 
     2 (D[r[u, v], u] D[r[u, v], v])/r[u, v]^2 + 
     1/2 D[F[u, v], u, v]/F[u, v] + D[r[u, v], u, v]/r[u, v]);
sol = ParametricNDSolve[{D[f[u, v], u, 
      v] == -(W[u, v, 2]) f[u, v], 
    D[g[u, v], u, v] == -(W[u, v, 2]) g[u, v], f[u, -0.001] == 10, D[f[0.001, v], v] == dfdu, 
    g[u, -0.001] == -8, D[g[0.001, v], v] == dgdu}, {f, g}, {u, 0.001,
     200}, {v, -200, -0.001}, {dfdu, dgdu}, PrecisionGoal -> 0, 
   WorkingPrecision -> MachinePrecision(,MaxStepSize[Rule]1)]
  dfsol = FindMinimum[((f[dfdu, dgdu] /. sol)[200, -0.001])^2 // 
  Evaluate, {{dfdu, 1}, {dgdu, 1}}]
  Plot3D[(f[dfdu, dgdu] /. sol /. dfsol[[2]])[u, v] // Evaluate, {u, 
   0.001, 120}, {v, -120.1, -0.001}, PlotRange -> All]

I get a very long list of errors saying that Solve uses inverse functions, that initial conditions are inconsistent and more.

The problem (I think) is that I am not specifying the initial conditions correctly I think I sould use $\dfrac{df}{du}\Bigg|_{0.001,-0.001}$ and $\dfrac{dg}{du}\Bigg|_{0.001,-0.001}$, but I do tno know how to write it properly.

Can anyone help?

EDIT I have solved part of the problem: the condition on the derivative should be D[f[u,-0.001],u]=dfdu and analogously for the other, but now I get the error ``

ParametricNDSolve::fembdnl: The dependent variable in f$4752==1. in > the Dirichlet boundary condition DirichletCondition[f$4752==1.,v$4727==-0.001] needs to be linear. >>