function F = coulombforce(q1,q2,x1,y1,z1,x2,y2,z2)
    k = 9*10^9;
    R = sqrt((x2-x1)^2+(y2-y1)^2+(z2-z1)^2);
    F = k*q1*q2/R^2;
end