 
 
 //  ~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
-            Fac  = fmin(FacMax,fmax(FacMin,FacSafe/pow(Err,ONE/ros_ELO)));
+            Fac  = fmin(FacMax,fmax(FacMin,FacSafe/pow(Err,ONE/4.0)));
             Hnew = H*Fac;
 
 //  ~~~>  Check the error magnitude and adjust step size



__device__  static  int ros_Integrator_ros4(double * __restrict__ var, const double * __restrict__ fix, const double Tstart, const double Tend, double &T,
        //  Rosenbrock method coefficients
        const double * __restrict__ ros_M, const double * __restrict__ ros_E, const double * __restrict__ ros_A, const double * __restrict__  ros_C, 
        const double * __restrict__ ros_Alpha, const double * __restrict__ ros_Gamma, const double ros_ELO, const int * ros_NewF, 
        //  Integration parameters
        const int autonomous, const int vectorTol, const int Max_no_steps, 
        const double roundoff, const double Hmin, const double Hmax, const double Hstart, double &Hexit, 

            {
                int istage=1;
                for (int i=0; i<NVAR; i++){		
                    varNew(index,i) = K(index,0,i)*(2.0)  + var(index,i);
                }

                Fun(varNew, fix, rconst, varNew, Nfun,VL_GLO); // FCN <- varNew / not overlap 

                HC = (-0.7137615036412310E+01)/(direction*H);
                for (int i=0; i<NVAR; i++){
                    K(index,1,i) = K(index,0,i)*HC + varNew(index,i);
                }

                if ((!autonomous))
                {
                    HG = direction*H*(-0.1769193891319233E+01);
                    for (int i=0; i<NVAR; i++){
                        K(index,istage,i) += dFdT(index,i)*HG;
		     }
