I've already solved the problem. In line 351 I called "aceleracionx" twice instead of using the y component, now I get the circular orbits perfectly! :D
I'm supposed to give them random velocities as well. But, just for the sake of checking, I'm giving them the orbital speed (not random at all), and I do not get a circular orbit. Therefore, there has to be a problem.
If you read the code, line 87, there's the formula
$$v=\sqrt{\frac{GM}{r}}$$...
To simulate the trayectories of solar systems around a black hole (i.e. a galaxy) I have 3 classes in C++: cSystem, cBlackHole and cGalaxy. cSystem assigns initial values of position, velocity, etc to a solar system. cBlackHole does the same but just for the black hole. And cGalaxy mixes...
Hey, I'm a physics undergraduate (currently in 4th year) in Spain. I felt like these forums are great for discussion, and also for help. I'm really happy to be a part of this community. Hope to see y'all again!