You need to rethink your problem statement and initial conditions. The initial positions can't possibly be in the neighborhood of 100 meters since that is practically at the center of the Earth. And initial velocities in the neighborhood of 0.1 m/s don't make sense either because that would effectively produce nearly rectilinear motion with everything falling almost straight down. Are these supposed to be "relative" positions and velocities about some reference circular orbit at some altitude?

Also not sure why you are running an N-Body problem in the first place instead of N separate 2-Body problems. Are you trying to simulate the effects of the weak gravitational interaction between 4 small bodies as they orbit in a close group above a much larger Earth? Or ...?

You list the mass of each small body as 10^10 kg, but if I did my calculations correctly a spherical rock with that mass would have a radius of about 98 meters. That means these small bodies are practically touching each other. Is that the intent?

Your current description with the initial positions, velocities, and masses simply does't make sense and you need to correct this and post an update to your question.

SIDE NOTE: It is hard to follow your indexing, but at first glance the velocity derivatives don't look correct. I would have expected to see cumulative vector additions to calculate these derivatives, but I don't see that. I might suggest that you reshape your state vector upon entry to be 3x2N matrix to make your indexing easier. E.g.,

Then the first body position is simply Y(:,1) and its velocity is Y(:,2). The second body Pos & Vel would be Y(:,3) and Y(:,4). The n'th body Pos & Vel is Y(:,2*n-1) and Y(:,2*n). Etc. This will be much easier to write code downstream than what you currently have.

Alternatively you could reshape to a 6xN matrix, e.g.,

Then the first body position is Y(1:3,1) and the velocity is Y(4:6,1). The second body Pos & Vel is Y(1:3,2) and Y(4:6,2). The n'th body Pos & Vel is Y(1:3,n) and Y(4:6,n). Etc.

I suppose you could even get cute and create two helper indexing variables:

Then the n'th body position and velocity would simply be Y(pos,n) and Y(vel,n) which would make your downstream code even more readable.

E.g., using this last scheme you can get all the position derivatives in one line:

Then run your loops to fill in the dYdt(vel,:) parts.

After calculating the 3x2N or 6xN matrix derivative dYdt, then reshape it to a column vector for return purposes: