A simple orbital model might be developed like this:
day = 86400; G=6.67408*10^-11; ER=6367440; earthmass=5.9722*10^24; pot[{x_, y_, z_}] := -G earthmass / Sqrt[x^2 + y^2 + z^2]; grad = -Grad[pot[{x[t], y[t], z[t]}], {x[t], y[t], z[t]}]; odesys = { {x[0], y[0], z[0]} == {1.1 ER, 0, 0}, {x'[0], y'[0], z'[0]} == RotationTransform[ 98 \[Degree], {1, 0, 0}]@(7543.7 {0, 1, 0}), {x''[t], y''[t], z''[t]} == grad} ; soln = NDSolve[odesys, {x, y, z}, {t, 0, 5 day} ]; And that works just fine:
Show[ { Graphics3D[ {Opacity[0.5], Sphere[{0, 0, 0}, ER]}], ParametricPlot3D[{x[t], y[t], z[t]} /. soln, {t, 0, 5 day}, PlotRange -> Table[{-2 ER, 2 ER}, 3], AspectRatio -> Automatic] } ] 
But what I would like to do is just use the variables as vectors, something like:
odesys = {x[0] == {1.1 ER, 0, 0}, x'[0] == {0, 7534, 0}, x''[t] == -earthmass Norm[x[t]]^-2 G Normalize[x[t]]} which works
NDSolve[odesys, x, {t, 0, 5 day}] But I'm not sure how to go about getting Grad to take a variable and make the assumption that it is a 3-element vector. How should I go about this?
