<<PhysicalConstants`
kappa=GravitationalConstant[[1]];
casTed=AbsoluteTime[Date[]]
Δt=120;
casMinus=casTed-Δt;
casPlus=casTed+Δt;

SetAttributes[symbol,Listable];
symbol[planeta_String]:=ToExpression[StringJoin["r",planeta]];

objekty={"Earth","Sun","Venus"};
ClearAll[pocPoloha];
pocPoloha[planeta_String]:=Module[{rr,poloha},
rr=StringJoin["r",planeta];
poloha=AstronomicalData[planeta,{"Position",DateList[casTed]}];
ToExpression[rr][0]==poloha
];
pocPoloha[planety_List]:=pocPoloha/@planety;
pocPodm1=pocPoloha[objekty];

ClearAll[pocRychlost];
pocRychlost[planeta_String]:=Module[{rr,rychlost},
rr=StringJoin["r",planeta];
rychlost=(1/(2*Δt))(AstronomicalData[planeta,{"Position",DateList[casPlus]}]-AstronomicalData[planeta,{"Position",DateList[casMinus]}]);
ToExpression[rr]'[0]==rychlost
];
SetAttributes[pocRychlost,Listable];
pocPodm2=pocRychlost[objekty];

hmotnost[obj1_String]:=hmotnost[obj1]=AstronomicalData[obj1,"Mass"];

ClearAll[sila];
sila[dvoj:{obj1_String,obj1_String}]:=0;
sila[dvoj:{obj1_String,obj2_String}]:=Module[{r1,r2,m1,m2},{r1,r2}=symbol[dvoj];
m1=hmotnost[obj1];
m2=hmotnost[obj2];
kappa*((m1*m2)/((Norm[-r1[t]+r2[t]])^3))*(-r1[t]+r2[t])
];

SetAttributes[rovnice,Listable];
rovnice[obj_String]:=Module[{m,r,sily},
m=hmotnost[obj];
r=symbol[obj];
sily=sila[{obj,#}]&/@objekty;
m*r''[t]==Total[sily]
]
rovs1=rovnice[objekty]

ClearAttributes[Total,Listable]
Total[{((7.925926899061922`*^44 (rEarth[t]-rSun[t]))/((Norm[rEarth[t]-rSun[t]])^3)),{0,0,0}}]

Total[sila[{"Earth",#}]&/@objekty]

planetySolve[tVyp_?Positive]:=Module[{rovs},
rovs=Flatten[{rovs1,pocPodm1,pocPodm2}];
NDSolve[rovs,symbol[objekty],{t,0,tVyp}]
]

tmax=4*3600*24*365.25;
res1=planetySolve[tmax][[1]]

symbol[objekty]

{rEarth,rSun,rVenus}

ParametricPlot3D[rEarth[t]/.res1,{t,0,tmax}]