function dYdt = kepler_ode_manyBodies(t,Y)
% assume 'Y' is a vector of form Y =[x1; y1; vx1; vy1; x2; y2; vx2; vy2; ... xn; yn; vxn; vyn]
%1 = Sun, 2 = planet, 3...n = asteroids
% use the following four variables from the main worksheet
% must declare them 'global' here and there
global ms mp ma nAsteroids
rs = [Y(1);Y(2)];
vs = [Y(3);Y(4)];
rp = [Y(5); Y(6)];
vp = [Y(7); Y(8)];
aPosMatrix = zeros(2, nAsteroids);
aVelMatrix = zeros(2, nAsteroids);
for i = 1:nAsteroids
    curx = 5+4*i;
    cury = 6+4*i;
    curvx = 7+4*i;
    curvy = 8+4*i;
    aPosMatrix(1, i) = Y(curx, 1);
    aPosMatrix(2, i) = Y(cury, 1);
    aVelMatrix(1, i) = Y(curvx, 1);
    aVelMatrix(2, i) = Y(curvy, 1);
end


% Sun
Fsp = force12(rs, rp, ms, mp);%- ms * mp * G / norm(rp - rs)^3 .* (rp-rs);
totalFsa = 0;
for i = 1:nAsteroids
    totalFsa = totalFsa + force12(rs, aPosMatrix(:, i), ms, ma);
end
as = (Fsp + totalFsa) ./ ms;
%as = (Fsp+Fsa) ./ ms;
Fps = force12(rp, rs, mp, ms);%-ms * mp * G/ norm(rs - rp)^3 .* (rs - rp);
%Fpa = force12(rp, ra, mp, ma);%-ms * mp * G/ norm(ra - rp)^3 .* (ra - rp);
totalFpa = 0;
for j = 1:nAsteroids
    totalFpa = totalFpa + force12(rp, aPosMatrix(:, j), mp, ma);
end
ap = (Fps + totalFpa) ./ mp;
aAccMatrix = zeros(2, nAsteroids);

for k = 1:nAsteroids
    curaPos = aPosMatrix(:, k);
    Fas = force12(curaPos, rs, ma, ms);
    Fap = force12(curaPos, rp, ma, mp);
    aAccMatrix(:,k)  = (Fas + Fap) ./ ma;
end
%aAccMatrix
%pause()
aPosMatrix = aVelMatrix;
aVelMatrix = aAccMatrix;
%size(newA)
%size(aPosMatrix)
A = zeros(4, nAsteroids);
A(1:2, :) = aPosMatrix;
A(3:4, :) = aVelMatrix;

%creating dYdt
newY = zeros(8+nAsteroids*4, 1);
newYsp = [vs(1); vs(2); as(1); as(2); vp(1); vp(2); ap(1); ap(2)];  % must be column vector 

newY(1:8, 1) = newYsp;
cur = 8;
for m = 1:nAsteroids
    for p = 1:4
        cur = cur + 1;
        newY(cur, 1) = A(p, m);
    end
end

dYdt = newY;
