First, some slight changes to `orbita`: orbita[a_?NumberQ, b_?NumberQ, c_?NumberQ, d_?NumberQ, e_?NumberQ] := (orbita[a, b, c, d, e] = {X[t], Y[t]} /. First@NDSolve[{Vx'[t] == (-e*G*X[t]*Msun)/(X[t]^2 + Y[t]^2)^(3/2), Vy'[t] == (-e*G*Y[t]*Msun)/(X[t]^2 + Y[t]^2)^(3/2), X'[t] == Vx[t], Y'[t] == Vy[t], Vx[0] == a, Vy[0] == b, X[0] == c, Y[0] == d}, {Vx, Vy, X, Y}, {t, 0, 10^4}]); Memoization (or caching; see [`FindFit`](http://reference.wolfram.com/mathematica/ref/FindFit.html) documentation: [1](http://reference.wolfram.com/language/ref/FindFit.html#246277593) [2](http://reference.wolfram.com/language/ref/FindFit.html#703708898)) is done with `=`; the `:=` causes the return value to be `Null` and is the source of your first error message. The other problem is the what to return. I suggest `{X[t], Y[t]}`, with the variable `t` in place. I'd also strip an extra set of braces `{}` with `First@NDSolve...`. Second, I don't think `FindFit` will work on 2D univariate data (`{x, y}` as a function of `t`). At least I could find no example and a naive toy trial failed. So use `FindMinimum` to minimize the sum of squares. The objective function is given by ClearAll[obj]; obj[Mc_?NumericQ] := Total[#.# & /@ ((orbita[0, 4*10^7, 1000*UA, 500*UA, Mc] /. t -> # & /@ tempi) - dati)] Then minimize: {min, sol} = FindMinimum[obj[Mc], {Mc, 10^6}] > FindMinimum::lstol: The line search decreased the step size to within the tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the function. You may need more than MachinePrecision digits of working precision to meet these tolerances. >> (* {1.54869*10^9, {Mc -> 4.3*10^6}} *) Inspect the solution to see if the warning is significant. Show[ ParametricPlot[ orbita[0, 4*10^7, 1000*UA, 500*UA, Mc] /. sol, {t, 0, Max[tempi]}, PlotStyle -> None, Mesh -> {tempi}, MeshStyle -> {PointSize[Large], Red}], ListPlot[dati]] ![Mathematica graphics](https://i.sstatic.net/8EwEl.png) It does not look too bad.