// n - число переменных 

void fcn( int n, double ymas[],double f[]);

void fcn( int n, double ymas[], double f[])
{
	double x,y,z,vx,vy,vz,r,r2,r3,r5;
	double mu=3.9860044e5, // км^3/c^2
	       w,RE=6378.136,  // км - экваториальный радиус Земли
	       R0,h0=35786,     // км - высота геостационара
	       J2 = 1082.628e-6; // коэффициент второй зональной гармоники
	       R0=RE+h0;        
	       w=sqrt(mu/(R0*R0*R0)); // c^-1 - угловая скорость движения по круговой орбите без возмущений
// введение удобных обозначений для проверки правильности написанного:
	  x=ymas[0];
	  y=ymas[1];
	  z=ymas[2];
	  vx=ymas[3];
	  vy=ymas[4];
	  vz=ymas[5];
	  r2=x*x+y*y+z*z;
	  r=sqrt(r2);
	  r3=r*r2;
	  r5=r3*r2;
// собственно вычисление правой части:
	  f[0]=vx;        
      f[1]=vy;
      f[2]=vz;
      f[3]=-mu*x/r3+(1.5*J2*mu*RE*RE/r5)*x*(5.0*z*z/r2-1.0)+w*w*x+2.0*w*vy;
      f[4]=-mu*y/r3+(1.5*J2*mu*RE*RE/r5)*y*(5.0*z*z/r2-1.0)+w*w*y-2.0*w*vx;
      f[5]=-mu*z/r3+(1.5*J2*mu*RE*RE/r5)*z*(5.0*z*z/r2-3.0);
  return ;
}