Q:

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.

