Numerous bodies orbiting loop help

I've recently finished making a simulation of one body orbiting a stationary one. However I would like to extend it to more than 2 bodies like 50 or more. I would like to use a for loop instead of copying bits of code 50 times over. However this is where I'm stuck. I have an idea of what to do but I don't know how to code it.Any help will be appreciated :)

This is my function code:

function [ de ] = orbit(t,y)
de = zeros(4,1);
%Position
xx=y(1);
yy=y(2);
%Radius
r=(xx.^2+yy.^2).^0.5;
%Constants
M = 4E30; % mass
G = 6.67384E-11; 
%dX/dt
de(1) = y(3); %vx
de(3) = (-G.*M.*xx)/(r.^3); %ax
%dY/dt
de(2) = y(4); %vy
de(4) = -G.*M.*yy/(r.^3); %ay
end

and call code:

x0=149.513e9;
y0=0;
vx0=0;
vy0=29.78e3;
trange=[0:1:100]; %time range to be solved for
p0=[x0;y0;vx0;vy0];  %assemble an intial p
[tarray, parray] = ode45(@orbit,trange,p0);
x=parray(:,1);
y=parray(:,2);
plot(x,y)
xlabel('x');
ylabel('y');
title('y=f(x)');

[James Tursa]

Do you want to store all of the orbit results in one big matrix, or are you only interested in plotting the results (e.g., using hold on) as they are generated?

[sweetdreams]

I would like to see the numerical results when possible but the main thing is to be able to see all the bodies orbiting

[Walter Roberson]

Your posting is incomplete for the same reason that your previous exactly-the-same postings were incomplete: you do not indicate what you want to be looped over. For example should exactly the same equations be used each time but you want to use different initial positions? Or do the masses change?

[sweetdreams]

so what I was thinking was that I would define something like 'noplanets = 48' in the call code with 4x noplanet initial condition [x,y,vx,vy] then stack them for however number of sets of 4 there are [x1,y1,vx1,vy1,x2,y2...]. Then in the orbit code, split y into sets of 4 and for each set there are 4 derivatives and stack the output. And find noplanet size(y)/4. E.g. for loop 1 it will be elements 1->4 so y(1:4). for loop 2 it will be elements 5->8 so y(5:8) so for ith element it will be y(4i-3:4i) and take out set of 4 derivatives to be called and plotted. mass will not change as it is the mass of the 'black hole' which is at the origin. the only thing that will change will be the initial conditions.

Answer 1

[Walter Roberson]

function [ DE ] = orbit(t,Y)
  y = reshape(Y,4,[]);
  de = 0 * Y;
  %Position
  xx=y(1,:);
  yy=y(2,:);
  %Radius
  r=(xx.^2+yy.^2).^0.5;
  %Constants
  M = 4E30; % mass
  G = 6.67384E-11; 
  %dX/dt
  de(1,:) = y(3,:); %vx
  de(3,:) = (-G.*M.*xx)./(r.^3); %ax
  %dY/dt
  de(2,:) = y(4,:); %vy
  de(4,:) = -G.*M.*yy./(r.^3); %ay
  DE = de(:);
end