Determing Velocity of Moons

I have a question that I believe is relatively easy to answer, I am working on an $N$-body simulation of a fictional star system and am having trouble finding the velocity of moons so that they will remain in near-circular orbits (were they not to be influenced by other moons). So far, I can calculate the velocities of planets relative to their parent star in order to have a near circular orbit at a specific semi-major axis. When I say velocity, I’m looking for a 3-dimensional vector and not just a 1D value. I find the velocity of the planets by calculating the orbital speed like$$
v_0 ~approx~ sqrt{frac{GM}{r}}
,,$$where $r$ is the length of the planet’s semi-major axis and $M$ is the mass of the star. This value is then multiplied by a 3-dimensional vector that points tangential to the planet’s desired orbit which gives me the velocity. This works fine but I’m having trouble doing this for moons of these planets. I made an attempt by repeating the process that I used for the planets but instead $r$ becomes the moon’s semi-major axis from its parent planet and $M$ becomes the mass of the planet. I then find the orbital velocity and multiply it by a vector tangential to the desired orbit relative to the planet and add the planets orbital velocity (as a 3D vector) which doesn’t work as many moons get left behind, others fall into the planet, and few remain in stable orbits (without a lot of perturbations from other objects).

While the moon’s velocity relative to the planet (so that it would have a circular orbit) is being properly calculated, the moon’s velocity relative to the star isn’t as apparently simply adding the planet’s velocity doesn’t work – and this is what I’m looking for. Can anyone help me with finding the velocity (relative to the sun) of a moon so that it will remain in a stable orbit around its parent planet?

EDIT: Code that calculates velocity of moons

double localOrbitalSpeed = sqrt((g * parentObjectMass)/(localSMA));
velocity = normalize(cross(position, Vector3(0, -1, 0))) * localOrbitalSpeed + parentObjectVelocity; 

and the code for planets

double orbitalSpeed = sqrt((g * starMass)/(starSMA));
velocity = normalize(cross(position, Vector3(0, -1, 0))) * orbitalSpeed;

and that’s pretty much all that really determines initial velocities, not much room for… while writing this I figured out what I was doing wrong, when calculating the velocity vector, the tangent vector is being determined relative to the star instead of the planet :/ so yeah, simple programming error. I just fixed this and will see if I have any more issues.

EDIT 2: Well now that the velocity is being properly calculated, I find that the moons tend to spiral into the planet or spiral away in only a few orbits (still without influencing each other), I know that moons do, do this but due to tidal forces (which I am not simulating), right? I assume that the integration method I used is causing this problem.