You could probably expect that $\langle v_r\rangle=\langle v_z\rangle=0$--that is, most of your velocities will be in the $\hat\phi$ direction (assuming you are using cylindrical coordinates).1
For stable orbits, the kinetic and potential energies must be equal, so you should end up seeing that
$$
v(r)=\sqrt{-2\Phi(r)}
$$
where $\Phi$ is the gravitational potential energy (and there are a few choices). You also can add Normally-distributed values to the velocity components (as suggested here) to give a little perturbations to the orbits. It is pretty straight-forward to convert to another coordinate system from here, if need be.
I also discuss in this related question an algorithm to pick the velocities assuming a Plummer model (discussed in the first link in the post). I expect a similar algorithm could be developed for the alternative potentials, but haven't worked the math of it.
Note also that, while the system of equations for the $N$-body problem is pretty straight-forward,
\begin{eqnarray}
\mathbf v&=\frac{\mathrm d\mathbf x}{\mathrm d t}\\
\mathbf F&=m\frac{\mathrm d\mathbf v}{\mathrm dt}
\end{eqnarray}
writing it in code is actually quite hard for large $N$. The force $\mathbf F$ is computed between each pair of objects in the $N$-body system:
$$
\mathbf F_i=-\sum_{j\neq i}G\frac{\hat{\mathbf x}_{ij}}{\left|\mathbf x_{ij}\right|^2}
$$
where $\mathbf x_{ij}=\mathbf x_i-\mathbf x_j$ and $\hat{\mathbf{x}}_{ij}$ is the direction of the force between bodies $i$ and $j$. This means you automatically have at least an $\mathcal O(N^2)$ problem (at least in time, it's $\mathcal O(N)$ in memory).2
A real galaxy has $N\sim10^{11}$ stars, which is probably not possible for any modern computer to handle.3 I believe that $N\sim10^6$ is a good round number for modelling galaxies, but these can still take hours on computer clusters (depending on what the simulation does) and I'll go out on a limb and guess that you don't have a computer cluster available to you, so you may want to try $N\sim10^3$ or $N\sim10^4$ instead--but don't be surprised if this still takes a very long time!
1. It is probably easiest to use $r,\phi,z$ coordinates for generating the velocities, then transforming them back to Cartesian (assuming you are working in Cartesian).
2. There are tricky algorithms that can reduce this to $\mathcal O\left(n\log n\right)$, but this may be a bit much for your purposes.
3. Though there was the one trillion body simulation, they used the whole K-computer for something like a month straight (if I recall the news articles correctly). I doubt you have that type of resources available.
