The $N$-body problem is an initial value problem, i.e., given a system in a given state at time $t=t_0$, what will its state be at an arbitrary time $t$?
For the $N$-body (gravitational) problem, the question would be: given a system with $N$ masses, each with a given position and velocity, where will they be at a time $t$ in the future?
Thus, to be able to determine (numerically) the positions of $N$ bodies interacting under gravity, their initial positions and velocities must be known.
Given these initial positions and velocities, you can then start to time-step the simulation forward using an integrator of your choosing.
As you have stated, the force between two masses interacting under gravity is determined by the distance between them. This is given by Newton's law of gravitation:
$$F=\frac{Gm_1m_2}{r^2}\tag{1}$$
In vector form, the force on body 1 due to body 2 is given by:
$$\textbf{F}_{12}=-\frac{Gm_1m_2}{r_{12}^2}\hat{\textbf{r}}_{21}\tag{2}$$
... where $\hat{\textbf{r}}_{21}$ is the unit vector pointing from body 2 to body 1.
This is trivially generalised to $N$ bodies, where the force on a body $i$ due to all other bodies is:
$$\textbf{F}_i=-Gm_i\sum_{\substack{j=1 \\ j \neq i}}^N{\frac{m_j}{r_{ij}^2}\hat{\textbf{r}}_{ji}}\tag{3}$$
... which can be divided by $m_i$ to give the acceleration experienced by body $i$:
$$\textbf{a}_i=-G\sum_{\substack{j=1 \\ j \neq i}}^N{\frac{m_j}{r_{ij}^2}\hat{\textbf{r}}_{ji}}\tag{4}$$
Given all of the initial positions and velocites of the bodies, the initial acceleration $\textbf{a}_i$ for each body will be known (exactly), and can then be used to approximate the velocity one time-step later. Conversely, the initial velocity of each body (which is also known) can be used to approximate the position one time-step later.
Using the Euler method of integration (probably the worst, but also the simplest to grasp), the time-stepping would look like this:
$$\textbf{r}_{n+1} = \textbf{r}_n + \textbf{v}_n\Delta{t}\tag{5}$$
$$\textbf{v}_{n+1} = \textbf{v}_n + \textbf{a}_n\Delta{t}\tag{6}$$
... where $\textbf{r}_n$, $\textbf{v}_n$ and $\textbf{a}_n$ are the position, velocity and acceleration vectors of a body at the start of the time-step, respectively, $\textbf{r}_{n+1}$ and $\textbf{v}_{n+1}$ are the position and velocity vectors at the end of the time-step, respectively, and $\Delta{t}$ is the size of the time-step itself. $\textbf{a}_n$ is determined by plugging all the positions of the bodies at the start of the step into $(4)$.
The global error per step for the Euler method is $\mathcal{O}(\Delta{t})$ (not so good).
Thus, to elaborate on your question, in the case of the Euler method, you would simply plug your initial positions and velocities of all the particles into equations $(5)$ and $(6)$ (having already calculated their accelerations $\textbf{a}_i$ using $(4)$) to obtain approximations for their positions and velocities at a time $\Delta{t}$, and then you would repeat the procedure for as many iterations as you like to obtain the positions and velocities of the particles at a time $2\Delta{t}$, $3\Delta{t}$, $4\Delta{t}$, ... so on and so forth.
I hope this helped.
