El problema de Euler de los tres cuerpos

Sistema de ecuaciones diferenciales de segundo orden

Resolver el sistema de ecuaciones diferenciales de segundo orden

$\begin{array}{l} \frac{d^{2} x}{d t^{2}} = - 4 π^{2} (\frac{x}{{(x^{2} + y^{2})}^{3 / 2}} + \frac{α (x - d)}{{({(x - d)}^{2} + y^{2})}^{3 / 2}}) \\ \frac{d^{2} y}{d t^{2}} = - 4 π^{2} (\frac{y}{{(x^{2} + y^{2})}^{3 / 2}} + \frac{α y}{{({(x - d)}^{2} + y^{2})}^{3 / 2}}) \end{array}$

Definimos una nueva energía e por unidad de masa en este sistema de unidades

$e = \frac{1}{2} (v_{x}^{2} + v_{y}^{2}) - 4 π^{2} (\frac{1}{\sqrt{x^{2} + y^{2}}} + \frac{α}{\sqrt{{(x - d)}^{2} + y^{2}}})$

Con las condiciones iniciales siguientes: en el instante t=0, x=x₀, y=y₀, v_x=v_0x, v_y=v_0y.

public class Sistema extends RungeKutta{
    double alfa; 
    double d;
    public Sistema(double alfa, double d,  double h){
      super(h);
      this.alfa=alfa;
      this.d=d;
    }
    public double f(double x, double y, double vx, double vy, double t){
         double r1=Math.sqrt(x*x+y*y);
         double r2=Math.sqrt((x-d)*(x-d)+y*y);
         double z=-4*Math.PI*Math.PI*(x/(r1*r1*r1)+alfa*(x-d)/(r2*r2*r2));
         return z;
    }
    public double g(double x, double y, double vx, double vy, double t){
         double r1=Math.sqrt(x*x+y*y);
         double r2=Math.sqrt((x-d)*(x-d)+y*y);
         double z=-4*Math.PI*Math.PI*y*(1.0/(r1*r1*r1)+alfa/(r2*r2*r2));
         return z;
    }
    public double energia(double x, double y, double vx, double vy){
         double r1=Math.sqrt(x*x+y*y);
         double r2=Math.sqrt((x-d)*(x-d)+y*y);
       double z=(vx*vx+vy*vy)/2-4*Math.PI*Math.PI*(1.0/r1+alfa/r2);
       return z;
    }
}
//Objeto de la clase Sistema
void setNuevo(double alfa, double d, double x0, double y0, double V0x, double V0y){
    this.d=d;
    estado=new Estado(0.0, x0, y0, V0x, V0y);
    sis=new Sistema(alfa, d, calculaPaso());
   eInicial=oscilador.energia(x0, y0, V0x, V0y);
     error=0.0;
 }
 double calculaPaso(){
     double r1=Math.sqrt(estado.x*estado.x+estado.y*estado.y);
     double r2=Math.sqrt((estado.x-d)*(estado.x-d)+estado.y*estado.y);
     double distMed=Math.sqrt(r1*r2);
     double dt=0.0001*Math.pow((distMed/(d*0.1)), 1.5);
     return dt;
 }
 void mover(){
     sis.resolver(estado);
     sis.setPaso(calculaPaso());
    double energia=oscilador.energia(estado.x, estado.y, estado.vx, estado.vy);
   error=Math.abs((energia-eInicial)*100/eInicial);
   if(error>1.0){
      //se para;
   }
 }