void Function(double x, double *y, double *f)
{   double r12, r13, r23;
    
    r12 = (y[0]-y[3])*(y[0]-y[3])+(y[1]-y[4])*(y[1]-y[4])+(y[2]-y[5])*(y[2]-y[5]);
    r12 = pow(r12, 1.5);
    r13 = (y[0]-y[6])*(y[0]-y[6])+(y[1]-y[7])*(y[1]-y[7])+(y[2]-y[8])*(y[2]-y[8]);
    r13 = pow(r13, 1.5);
    r23 = (y[3]-y[6])*(y[3]-y[6])+(y[4]-y[7])*(y[4]-y[7])+(y[5]-y[8])*(y[5]-y[8]);
    r23 = pow(r23, 1.5);
	f[0] = y[9];
	f[1] = y[10];
	f[2] = y[11];
	f[3] = y[12];
	f[4] = y[13];
	f[5] = y[14];
	f[6] = y[15];
	f[7] = y[16];
	f[8] = y[17];
	f[9] = 3*(y[3]-y[0])/r12 + 5*(y[6]-y[0])/r13;
	f[10] = 3*(y[4]-y[1])/r12 + 5*(y[7]-y[1])/r13;
	f[11] = 3*(y[5]-y[2])/r12 + 5*(y[8]-y[2])/r13;
	f[12] = -1*(y[3]-y[0])/r12 + 5*(y[6]-y[0])/r23;
	f[13] = -1*(y[4]-y[1])/r12 + 5*(y[7]-y[1])/r23;
	f[14] = -1*(y[5]-y[2])/r12 + 5*(y[8]-y[2])/r23;
	f[15] = -1*(y[6]-y[0])/r13 - 3*(y[6]-y[3])/r23;
	f[16] = -1*(y[7]-y[1])/r13 - 3*(y[7]-y[4])/r23;
	f[17] = -1*(y[8]-y[2])/r13 - 3*(y[8]-y[5])/r23;

	return;
} // fEval