Two body problem in processing

Learning Processing has been a fun (and frustrating) experience for me. Coding in Java has me realizing how much Python spoils me. I keep getting NullPointerExceptions in Java because of how I deal with objects that I haven't initialized yet. I like that I have to declare the type for every variable. This organizes my thinking, and forces me to be more deliberate with my code. Processing is a powerful tool for visualization, but it's not as easy to use for doing math as Numpy. In Numpy adding two arrays together is trivial. So is performing operations. In Java, at least with regular arrays, you have to iterate through the array. Ugh. 

#Python
x = np.random.random((10,))
y = np.sin(x)
h = np.random.random((10,)) + x

//Java/Processing
float[] x = new float[10];
float[] y = new float[10]; 
float[] h = new float[10];
for (int i=0; i < 10; i++){
x[i] = random(10);
y[i] = sin(x[i]);  
h[i] = random(10) + x[i];
}

Anyways, I set about solving the two body problem in Processing. Normally when solving the two body problem, you do a coordinate transform, which ultimately ends up turning the problem into a 1 body problem. I want to turn my code into a three or four body code, so I stuck with solving things in cartesian coordinates. You might be thinking that I'm crazy, but I don't mind the extra lines of code for the simplicity of my approach. Let's dive in. The first thing to do is to set up the Lagrangian for our system. The idea of the two body problem is solve for the equations of motion of two particles who are mutually gravitationally attracted to each other. The Lagrangian for the system is as follows. $$ L = T-V\\ L = \frac{m_1}{2}(\dot{x_1}^2 + \dot{y_1}^2) + \frac{m_2}{2}(\dot{x_2}^2 + \dot{y_2}^2) - \frac{Gm_1m_2}{\sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2}} $$ By applying the Euler-Lagrange equations we can arrive at the following equation of motion for the four coordinates in our system. $$ \ddot{x_1} = -Gm_2(x_1-x_2)[(x_1 - x_2)^2 + (y_1 - y_2)^2]^{-3/2}\\ \ddot{y_1} = -Gm_2(y_1-y_2)[(x_1 - x_2)^2 + (y_1 - y_2)^2]^{-3/2}\\ \ddot{x_2} = Gm_1(x_1-x_2)[(x_1 - x_2)^2 + (y_1 - y_2)^2]^{-3/2}\\ \ddot{y_2} = Gm_1(y_1-y_2)[(x_1 - x_2)^2 + (y_1 - y_2)^2]^{-3/2} $$ Now we have to solve a system of nonlinear second order differential equations. We can simplify this problem a lot by turning each of the second order equations into a first order equation. I've talked about this in previous posts, but I'll go over it again briefly. Take, for instance, the first coordinate $x_1$. Let me create new variables $x_1^a$ and $x_1^b$ defined as follows. $$ x_1^a = x_1\\ x_1^b = \dot{x_1}\\ \dot{x_1}^a = \dot{x_1} = x_1^b\\ \dot{x_1}^b = \ddot{x_1} $$ Now I can rewrite my second order equation for $x_1$ as a system of two first order differential equations in terms of $x_1^a$ and $x_1^b$. Note that I've gone ahead and introduced the 'a' and 'b' notation for the other coordinates. $$ \dot{x_1}^a = x_1^b\\ \dot{x_1}^b = -Gm_2(x_1^a-x_2^a)[(x_1^a - x_2^a)^2 + (y_1^a - y_2^a)^2]^{-3/2} $$ With these equations in hand I can use any first order ODE method to step my equations through time. I used Runge Kutta in my code. I coded all this up in Processing, and I made a .gif for your enjoyment! Check it out below.

You might notice some decaying of the orbit. This is not physical, rather it represents a loss of stability in my ODE solver. This is even more pronounced when using a simple method like forward Euler. Also notice that I make no assumptions about the position or mass of the particle at the center of the frame -- I just defined it as having about 100x the mass of the second particle.