Stable disk for $N$-body simulation I want to simulate a Galaxy of $N$ particles which I have to generate first. What I have done so fare does not lead to the result I want to see. I probably did an error and don't see it. Here comes what I did:
1) Generate uniformly distributed particles inside a unit circle with mass $m_i = 1/N$:
For every point I do:
   R = 1;
   radius = R*sqrt(rand(0,1));
   phi = rand(0,2*pi);
   m = 1/N;
   x = radius*cos(phi);
   y = radius*sin(phi);

Doing that I get the following distribution:

2) Generate velocities such that the circle does not explode or implode:
I want to generate something where $F_{\text{Gravity}} = F_{\text{Centrifugal}}$
Since $F_{\text{Centrifugal}} = \frac{mv^2}{r}$ and $F_{\text{Gravity}} = G\frac{M(r) m}{r^2}$ I can equate both equations and get for the velocity $v$
$M(r)$ is the mass inside radius $r$. I can write $M_{\text{tot}} = \pi R^2$ and $M=\pi r^2$ (density $\rho=1$). From those equations I get $M(r) = \frac{r^2}{R^2} M_{\text{tot}}$. I can now put it back to $v$ and get:
$v = \sqrt{\frac{Gr^2M_{\text{tot}}}{rR^2}}$
Since $M_{\text{tot}} = G = R = 1$ I get:
$v = \sqrt{r}$
So, for every point I do:
   velocity = sqrt(radius);
   v_x = -velocity*sin(phi);
   v_y = velocity*cos(phi);

Then I get the following result:




Ok, it seems that the circle explodes somehow. Particles from the center move outwards. Is my calculation for $v$ wrong, or does it simply not work for my purpose? Do I need another velocity distribution? In the end, for large $t$, I can't see any spiral arms forming or any other structure. I tried $N=4000$ and $N=50000$. Do I need more particles? Did I something wrong?
 A: Your logic is correct, including the choice of tangential velocity as
 x = r*cos(phi);
 y = r*sin(phi);
 vx=-v*sin(phi);
 vy= v*cos(phi);

Where you (tried to) set v in centrifugal balance. However, the gravitational attraction of a disc is different from that of a sphere, in particular your assumption that the acceleration equals GM(r)/r^2 is wrong. The best way to find the correct form is to obtain the forces from your N-body model and take the azimuthal average of the radial acceleration in radial bins. Your initial model should satisfy the virial theorem, I.e. the kinetic energy should be -0.5 times the gravitational energy.
But how did you compute the forces and how did you integrate the equations of motion? You should avoid the exact N-body forces, because these diverge at close encounters and instead use softened forces. You essentially add for each particle-particle interaction a constant vertical offset h, the softening length. Set this to about 0.1 times the radius.
The time integration is best done using the leapfrog integrator. For example
for(I:particles)
    I.vel += 0.5*tau*I.acc
    I.pos += tau*I.vel
time += tau
compute_accelerations(particles);
for(I:particles)
    I.vel += 0.5*tau*I.acc

to integrate one time step of length tau. Of course, you must compute the initial accelerations before the first step. The step size tau should be small enough for the total energy to be conserved by one part in 1000 or better.
