In my physics class, I had to calculate the trajectory of a projectile that was fired (very fast) with $v_0$ in an angle off a planet (radius $R$, mass $M$) from the surface. The projectile would reach a maximum distance and start to circle around the planet. The question was whether it would hit the planet. We solved that analytically, but I wonder how to simulate such a problem numerically.So basically I start off with $$\vec{r}_0 = \begin{pmatrix} 0 \\ 0 \\ R \end{pmatrix}$$and$$\vec{v}_0 = v_0 \begin{pmatrix} \cos(30°) \\ 0 \\ \sin(30°) \end{pmatrix}$$The acceleration on the projectile is just$$ \vec{a}(\vec{r}) = \frac{\vec{F}(\vec{r})}{m} = G \frac{M}{r^3} \vec{r} $$I could implement a simple C++/Python program which basically does:$$ \vec{r}_{t+1} = \vec{r}_{t} + \vec{v}_{t} \cdot \mathrm dt $$$$ \vec{v}_{t+1} = \vec{v}_{t} + \vec{a}_{t} \cdot \mathrm dt $$But this seems a little unsuited as in either language I do not have natural vector classes.Would Octave or Mathematica be better suited for this task? How would I implement this problem to obtain a list of x, y, z coordinates in order to plot the trajectory in gnuplot?
Implementing Euler's method for initial value ODEs
software;ode
Just use Python + NumPy. Here is the codefrom math import pi, sin, cosfrom numpy import array, vstackfrom numpy.linalg import normfrom pylab import plot, savefigR = 10.M = 1.G = 1.phi = 30 * 180. / piv0_magnitude = 0.1r = array([0, 0, R])v = v0_magnitude * array([cos(pi), 0, sin(pi)])trajectory = []dt = 0.1t = 0while t < 100:    trajectory.append(r.copy())    a = -G * M * r / norm(r)**3    v += a * dt    r += v * dt    t += dttrajectory = vstack(trajectory)plot(trajectory[:, 0], trajectory[:, 2])savefig(trajectory.png)And here is the graph:I have also uploaded the above code and the graph at this gist. 