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,
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)?
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?