Trouble Making 3rd-Order Sympletic Integrator for Planitary N-Body Problem (A Hamiltonian System)

I am doing a solar-system simulation. I am using Ruth’s 3rd order sympletic integrator to avoid the problem of Energy Drift (which I had with RK4), but the the planets quickly leave orbit, and energy is by no means conserved (just like with RK4).

I applied this to the N-body problem with the following:

(KE=1/2mv^2)

I have implemented this into Fortran 2008, where x, a, v, p, and m are all vectors of length 30, which hold the x,y,z position, x,y,z acceleration, x,y,z velocity, x,y,z momentum, and m,m,m respectively for 10 separate bodies in the solar system (Planets + Sun + Pluto).

Acceleration on each body is calculated as the sum of a=GM/(r^2) for x,y,z for each other body on each body.

Here is the integration part of the code:

!----------Looping Through Time-----------
do while(t<365.250000d0) ! Length of simulation in days
    !----------Calculating Values-----------
    call calc_acc(masses,x,a)
    p1=p+(7.0d0/24.0d0)*h*m*a
    x1=x+(2.0d0/3.0d0)*h*p1/m
    call calc_acc(masses,x1,a)
    p2=p1+(3.0d0/4.0d0)*h*m*a
    x2=x1-(2.0d0/3.0d0)*h*p2/m
    call calc_acc(masses,x2,a)
    p=p2-(1.0d0/24.0d0)*h*m*a
    x=x2+h*p/m
    v=p/m
    t=t+h
    !----------Saving Values-----------
    do bodnum=1,10,1
        write((100+bodnum),*) t, x((1+3*(bodnum-1)):(3+3*(bodnum-1))), v((1+3*(bodnum-1)):(3+3*(bodnum-1)))
        write((200+bodnum),*) x((1+3*(bodnum-1))), x((2+3*(bodnum-1))), x((3+3*(bodnum-1)))
    end do 
end do 

The full program can be found here.

Please tell me what I am doing wrong.

to_add=(rj-ri)*big_g*masses(i)*masses(j)/((abs(rj-ri))**3.0d0)