Hi all, I am implementing the LM method. The problem is that the model function Y(x) is the solution of an ODE like
F(x,Y(p1,p2,p3,...),dY(p1,p2,p3,...)/dx,p1,p2,p3,...)=0 (p=parameters)
which doesn't have a closed solution, i.e., there's no way to find y=Y(x). Thus, all I can do is to evaluate the dependent variable "y" for a known value "x" and a fixed set of parameters, with a numeric integration algorithm...
BUT the LM algorithm needs information like dY/dp1, dY/dp2, etc...

Then my question is: how can I estimate the partial derivatives only with values [xo,Y(xo)]???

I've seen the formula:

dY/dp~ [Y(xo,...p1+t...)-Y(xo,...p1...)]/t, with t~0, or

dY/dp~ [Y(xo,...p1+t...)-Y(xo,...p1-t...)]/(2t), with t~0

... but these formulas need a very high precision of the function evaluations, which are numerical estimates and already have the intrinsic error of the integration algorithm!!

Any advice?
Should I try a completely different fitting algorithm (ie, LM definitely does not work without good estimates of the Hessian matrix)?

Thanks.