 		double Fac = 0.5; //Step reduction if too many iterations
 		int NewtonIter = 0;
 		double Theta = 0;
+#ifdef NEWTON_UNROLL
+		#pragma unroll NEWTON_UNROLL
+#endif
 		for (; NewtonIter < NewtonMaxit; NewtonIter++) {
 			RK_PrepareRHS(t, pr, H, y, F0, Z1, Z2, Z3, DZ1, DZ2, DZ3);
 			RK_Solve(H, E1, E2, DZ1, DZ2, DZ3, ipiv1, ipiv2);

            }

            NewtonDone = (NewtonRate * NewtonIncrement <= NewtonTol);
            if (NewtonDone) break;
            if (NewtonIter == NewtonMaxit - 1) {
            	//todo implement return codes
				y[0] = logf(-1);
				return;
            }
		}
		if (!NewtonDone) {
			H = Fac * H;

		NewtonDone = false;
        Fac = 0.5; // ! Step reduction factor if too many iterations
        double NewtonIncrement = 0;
#ifdef NEWTON_UNROLL
		#pragma unroll NEWTON_UNROLL
#endif
        for (NewtonIter = 0; NewtonIter < NewtonMaxit; NewtonIter++) {
        	//!~~~>   Prepare the loop-dependent part of the right-hand side
        	WADD(y, Z4, TMP);
