How to solve $\dot{\mathbf{v}} = - \frac{GM}{r^3}\mathbf{r}$ using differential eq-n solver in MATLAB I need to solve second order system with boundary conditions with help of differential equation solver. The system is $$\dot{\mathbf{v}} = - \frac{GM}{r^3}\mathbf{r}.$$
I tried to use MATLAB dsolve command like this:
syms x(t) y(t)   
eqns = [diff(x,t,2) == -1*x/sqrt(x*x+y*y)^3, diff(y,t,2) == -1*y/sqrt(x*x+y*y)^3];   
Dx = diff(x,t);
Dy = diff(y,t);
cond = [x(0) == 10, y(0) == 0, Dx(0) == 0, Dy(0) == 1];   
dsolve(eqns,cond)
Warning: Explicit solution could not be found. 

EDIT:
I think i need to use ode45 instead, something like this:
f = @(t,y)[y(3); y(4); -y(1)/sqrt(y(1)*y(1)+y(2)*y(2))^3; -y(2)/sqrt(y(1)*y(1)+y(2)*y(2))^3];
[ts,ys]=ode45(f, [0, 100], [10; 0; 0; 0.1]);
plot(ys(:,1), ys(:,2));

And it kinda works:

Very good tutorial which helped me a lot: Using Matlab for Higher Order ODEs and Systems of ODEs
 A: You can use
function ydot = orbit(~, y)
  %ORBIT: differential equation for 2D orbit
  %  Returns the vector of the derivatives of y at time t.
  ydot = [y(3:4); -1/norm(y(1:2)) * y(1:2)];
end

When you call ode45 (with @orbit) you may want to set the MaxStep option (with odeset) to get better accuracy.  Of course, you can also write an anonymous function directly.  
EDIT: I just saw your edit.  The non-closed orbit is probably the result of not controlling the maximum integration step.
A: Another way is building a matrix $3 \times N$, and a digital differentiation filter
$$\frac{1}{\Delta t}\begin{bmatrix}1&-1\end{bmatrix}$$
Denote the matrix representation $\bf D$ of it's convolution on a vectorization of said image $\bf v$. To be strict, we should also use some interpolation filter to make sure we measure differential and function on the same points in time:
$$\frac{1}{2}\begin{bmatrix}1&1\end{bmatrix}$$
Let us call the matrix representation of this interpolation $\bf L$. Now let us solve a sequence of norm problems. Norm 2 (resulting in least-squares) should probably be good enough for just a simple test. $${\bf v_0} = \min_{\bf v}\left\{\left\|\frac{\bf Lv}{d^3}-\bf Av\right\| + {\bf C_i}\|\bf v - v_i\|\right\}$$
Where $d^3 = |{\bf v}|^3$ updated after each new solved one. $\bf C_i$ is a matrix encoding the certainty of boundary (initial value) condition and $\bf v_i$ the values of said initial/boundary values.
This equation system will have a set of equations "grabbing ahold" of neighbouring time points implicitly giving a very powerful processing.
I'll get back if I get a chance to try it.
