function f = funtb(x,y)
  par = (y(1)^2+y(2)^2)^(3/2);
  f1 = -y(1)/par;
  f2 = -y(2)/par;
 f = [f1 f2]';
 end