nbody finite diff aprox

I am writing a program in matlab that will simulate the n-body problem.

The simulation is in 3 dimensions.

I want something like this:

run a script and then it prompts the user for how many bodies to simulate, assigns random masses, random velocities, and random positions.

It gets its new position using finite difference approximation.

I know I must use f=-gmM/r^2 and f=mx'' for this.

But in all honesty I have no idea how to start this project off.... Can i get some tips? starting points? is my project doable?