##Equazione Differenziale
 #sdot(1) = vx
 #sdot(2) = vy
 #sdot(3) = ax
 #sdot(4) = ay
function sdot=f(stato,t)
  g = 9.8;
  x  = stato(1);
  y  = stato(2);
  vx = stato(3);
  vy = stato(4);
  
  ay = - g * (x*x)/(x*x + y * y) - (vx * vx + vy * vy) * y / (x * x + y * y);
  ax =   g * (x*y)/(x*x + y * y) - (vx * vx + vy * vy) * x / (x * x + y * y); 
  
  sdot(1) = vx;
  sdot(2) = vy;
  sdot(3) = ax;
  sdot(4) = ay;
endfunction
  
t = 0:0.01:10;
s0(1) = 5;
s0(2) = 0;
s0(3) = 0;
s0(4) = 0;

res = lsode("f", s0, t);
posizione_x = res(:, 1);
posizione_y = res(:, 2);


