function ydot = ode(y,x)
global nElectrons q m k nProtons
yRow=y';
r=yRow(1:3*nElectrons);
Fefinal=zeros(1,3*nElectrons);
Fpfinal=zeros(1,3*nElectrons);
for i=0:nElectrons-1
    for j=0:nElectrons-1
        if j~=i
            Fe(3*i+1:3*i+3)=-(r(3*j+1:3*j+3)-r(3*i+1:3*i+3)).*(k*q^2/((norm(r(3*j+1:3*j+3)-r(3*i+1:3*i+3)))^3));
            Fefinal(3*i+1:3*i+3)=Fefinal(3*i+1:3*i+3)+Fe(3*i+1:3*i+3); %sum electric forces for every electron from every electron
        end
    end
end
for i=0:nElectrons-1
    for j=0:nProtons-1
        Fp(3*i+1:3*i+3)=(x(3*j+1:3*j+3)-r(3*i+1:3*i+3)).*(k*q^2/((norm(x(3*j+1:3*j+3)-r(3*i+1:3*i+3)))^3));
        Fpfinal(3*i+1:3*i+3)=Fpfinal(3*i+1:3*i+3)+Fp(3*i+1:3*i+3);  %sum forces for every electron from every proton
    end
end
acceleration=(Fefinal+Fpfinal)./m;   %determine acceleration of particle (including direction) due to all forces on that particle
ydotRow=[yRow(3*nElectrons+1:6*nElectrons) acceleration];
ydot=ydotRow';
end

