jmenaTeles={"Earth","Jupiter","Sun"};

polohaTelesa[t_, planeta_String]:=AstronomicalData[planeta,{"Position",DateList[t]}]

generujData[jmenaTeles_]:=Module[{casTed},
(casTed=AbsoluteTime[];{jmeno->#,hmotnost->AstronomicalData[#,"Mass"],r0->polohaTelesa[casTed+100,#],v0->(polohaTelesa[casTed+100,#]-polohaTelesa[casTed-100,#])/200}&/@jmenaTeles)
]

dataTeles=generujData[jmenaTeles]

(*er[string_]:=ToExpression[StringJoin["r",string]];
er["Earth"]//FullForm*)

m=hmotnost/.dataTeles;

jmena=jmeno/.dataTeles;

pocR=r0/.dataTeles;

pocV=v0/.dataTeles;

<<PhysicalConstants`
κ=GravitationalConstant[[1]];
sila[naco_Integer,zdroj_Integer]:=κ*m[[naco]]*m[[zdroj]]*(r[naco][t]-r[zdroj][t])/Norm[r[naco][t]-r[zdroj][t]]^3

ClearAll[sila];
sila[naco_Integer,zdroj_Integer]:=
Module[{m1,m2,dr},
m1=m[[naco]];
m2=m[[zdroj]];
dr=r[naco][t]-r[zdroj][t];
κ*((-m1*m2*dr)/((Norm[dr])^3))
]
sila[i_Integer,i_Integer]=0;

pocetTeles:=Length[dataTeles];

ClearAll[silaNa];
silaNa[naco_Integer]:=Total[sila[naco,#]&/@Range[pocetTeles]]
silaNa[1]

ClearAll[silaNa];
silaNa[naco_Integer]:=Sum[sila[naco,i],{i,1,pocetTeles}]
silaNa[1]

soustava=dejRovnici/@Range[pocetTeles];

promene=r/@Range[pocetTeles]

jmena

tmax=365*24*3600;
resNahr=NDSolve[soustava,promene,{t,0,tmax}][[1]]

promene/.resNahr

ParametricPlot[(r[1][t]/.resNahr)[[{1,2}]],{t,0,tmax}]