Solución numérica de las ecuaciones del movimiento

Ecuación diferencial 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{d^{2} y}{d t^{2}} = - 4 π^{2} \frac{y}{{(x^{2} + y^{2})}^{3 / 2}} \end{array}$

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

public class Planeta extends RungeKutta{
Planeta(double h){
	super(h);
}
public double f(double x, double y, double vx, double vy, double t){
	double r=Math.sqrt(x*x+y*y);
	double z=-4*Math.PI*Math.PI*x/(r*r*r);
	return z;
}
public double g(double x, double y, double vx, double vy, double t){
	double r=Math.sqrt(x*x+y*y);
	double z=-4*Math.PI*Math.PI*y/(r*r*r);
	return z;
}
public double energia(double x, double y, double vx, double vy){
	double r=Math.sqrt(x*x+y*y);
	double z=(vx*vx+vy*vy)/2-4*Math.PI*Math.PI/r;
	return z;
}
}

//Objetos de la clase Planeta

	Estado estado=new Estado(0.0, x, 0.0, 0.0, vy); //estado inicial
	Planeta planeta=new Planeta(0.005);
 	eInicial=planeta.energia(estado.x, estado.y, estado.vx, estado.vy);

//rutina que calcula la trayectoria paso a paso
	planeta.resolver(valor);
	double energia=planeta.energia(estado.x, estado.y, estado.vx, estado.vy);
	double error=Math.abs((energia-eInicial)*100/eInicial);
	if(error>1.0){
		//detiene el movimiento
	}