How does NDSolve`ProcessEquations works?

Hello!

I have to compute geodesics of a 2-dimensions Riemannian manifold, for a number of points. The general form of the equation is (see the joined notebook):

NDSolve[EQ[u,v,t,A,B]]

where Eq is the Euler-Lagrange equation and A and B are the coordinates of the points (boundary conditions for the geodesic).
For instance,

{A, B} = {{0.045945904968576215, 0.19803051715809528}, {9.543534531634336, 1.6504203767297156}}

in the example posted.

Even with a small number of points, computations are very long, and I thought it was possible to save time by using combinations of (NDSolve\ProcessEquations, NDSolve\Iterate, NDSolve\ProcessSolutions) for changing initial condititions with a lower computation cost.
Clearly, it doesn’t work in my case… But is it merely impossible, or did I misused these commands (completely new for me)?

Thanks for your help,
Claude

Did you try ParametricNDSolve ?

No! I’m not experienced in differential equation solving, unfortunately… I’ll try in a minute.

Thank you very much.

I tried ParametricNDSolve and ParametricNDSolveValue, but Mathematica now returns expressions depending on the four parameters, which I cannot evaluate (even when three parameters are fixed).
What can I do?

Thanks a lot, Claude

After taking a quick look at your code, I suspect you’re not using the result of ParametricNDSolve correctly. The documentation example is

In[1]:= sol = 
 ParametricNDSolve[{y'[t] == a y[t], y[0] == 1}, y, {t, 0, 10}, {a}]

Out[1]= {y -> ParametricFunction[ <> ]}

In[2]:= y1 = y[1] /.  sol

Out[2]= InterpolatingFunction[{{0., 10.}}, <>]

and you are doing something different.

Sol1 = solGlob /. {b1, b2, e1, e2} -> Flatten[{A, BB}];

Yesss! The right manner is:

Sol1 = {Apply[u, Flatten[{A, B}]], Apply[v, Flatten[{A, B}]]} /. solGlob;

Thanks!