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 documentation: 1 2) 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]]
It does not look too bad.
