function z = ydot3(t,vector,mass)
%vector is 12-vector with position and velocity for all masses
%mass is 3-vector with M1,M2,and M3
m1=mass(1);g=1;mg1=m1*g;
m2=mass(2);mg2=m2*g;
m3=mass(3);mg3=m3*g;
px1=vector(1);py1=vector(3);pz1=vector(5);vx1=vector(2);vy1=vector(4);vz1=vector(6);
px2=vector(7);py2=vector(9);pz2=vector(11);vx2=vector(8);vy2=vector(10);vz2=vector(12);
px3=vector(13);py3=vector(15);pz3=vector(17);vx3=vector(14);vy3=vector(16);vz3=vector(18);

dist1_2=sqrt((px2-px1)^2+(py2-py1)^2 + (pz2-pz1)^2);   %Distance between M1 and M2
dist1_3=sqrt((px3-px1)^2+(py3-py1)^2+(pz3-pz1)^2);   %Distance between M1 and M3
dist2_3=sqrt((px3-px2)^2+(py3-py2)^2+(pz3-pz2)^2);   %Distance between M2 and M3

z=zeros(1,18);
z(1)=vx1;
z(2)=(mg2*(px2-px1))/(dist1_2^3) + (mg3*(px3-px1))/(dist1_3^3); % v-prime-x1
z(3)=vy1;
z(4)=(mg2*(py2-py1))/(dist1_2^3) + (mg3*(py3-py1))/(dist1_3^3); % v-prime-y1
z(5)=vz1;
z(6)=(mg2*(pz2-pz1))/(dist1_2^3) + (mg3*(pz3-pz1))/(dist1_3^3); % v-prime-z1

z(7)=vx2;
z(8)=(mg1*(px1-px2))/(dist1_2^3) + (mg3*(px3-px2))/(dist2_3^3); % v-prime-x2
z(9)=vy2;
z(10)=(mg1*(py1-py2))/(dist1_2^3) + (mg3*(py3-py2))/(dist2_3^3); % v-prime-y2
z(11)=vz2;
z(12)=(mg1*(pz1-pz2))/(dist1_2^3) + (mg3*(pz3-pz2))/(dist2_3^3); % v-prime-z2
z(13)=vx3;
z(14)=(mg1*(px1-px3))/(dist1_3^3) + (mg2*(px2-px3))/(dist2_3^3); % v-prime-x3
z(15)=vy3;
z(16)=(mg1*(py1-py3))/(dist1_3^3) + (mg2*(py2-py3))/(dist2_3^3); % v-prime-y3
z(17)=vz3;
z(18)=(mg1*(pz1-pz3))/(dist1_3^3) + (mg2*(pz2-pz3))/(dist2_3^3); % v-prime-z3