function dYdt = kepler_ode(t,y)
% assume 'Y' is a vector of form y = [ x ;  y ; vx ; vy ]

%result
%A = 1 1 0 0 1 1 0 0 ...
%B = 0 0 1 1 0 0 1 1 ...
%where dydt = (circshift(y,-2) * A) + (Fmatrix(condensed) * C * B) = 
% where A is the position of values that dont change ie velocities
% and B is the indices of values that do change


A = ones(size(y,1),1);
B = zeros(size(y,1),1);
for i=1 : size(A,1) 
    if rem(i+1,4) == 0 || rem(i,4) == 0 % every element that is at inidce 4n-1 or 4n
        A(i) = 0;                       
        B(i) = 1;
    end
end

global m N G 
Fij = zeros(N,N,2);
% since in force Fij = -Fji
% the matrix can be computed for Fij then added to its negative transpose
% Fmatrix = Fij + -1 .* Fij.'
% cutting n^2 computations to NlogN
%      1     2     3     .     .     .     N
% 1 |  0  | F12 | F13 | F1. | F1. | F1. | F1N |
% 2 |-F12 |  0  | F23 | F2. | F2. | F2. | F2N |
% 3 |-F13 |-F23 |  0  | F3. | F3. | F3. | F3N |
% . |-F1. |-F2. |-F3. |  0  | F.. | F.. | F.N |
% . |-F1. |-F2. |-F3. |-F.. |  0  | F.. |F.N  |
% . |-F1. |-F2. |-F3. |-F.. | F.. |  0  | F.N |
% N |-F1N |-F2N |-F3N |-F.N |-F.N |-F.N |  0  |
% this matrix can then be flattened to display all the total net force
% acting on each object.

for j = 1:N
    for i = 1:j
        if i == j
            Fij(i,j,:) =[0 0];
        else
            ri = [y(1+4*(i-1)) y(2+4*(i-1))];
            rj = [y(1+4*(j-1)) y(2+4*(j-1))];
            Fij(i,j,:) =  -G*m(i)*m(j)/(( norm(ri-rj) )^3).*   ( ri-rj );
        end
        
    end
end
Fmatrix = Fij + -1 .* (permute(Fij,[2 1 3]) );
%  Fmatrix(:,:,1) x
%  Fmatrix(:,:,2) y
Fmatrix = -1.*sum(Fmatrix);
%summation step/ flattening

amatrix = [Fmatrix(1,:,1) ./ (m'); Fmatrix(1,:,2) ./ (m')];
% |x|
% |y|
% inputing elements into be in the form of[0;0;ax1;ay1;0;0;ax2;ay2;...]
for i = 1:N
    B(3 + 4*(i-1)) = amatrix(1,i);
    B(4 + 4*(i-1)) = amatrix(2,i);
end

dYdt = circshift(y,-2).*A+B;









