 			if (NewtonIter > 0) 
 			{
 				Theta = NewtonIncrement / NewtonIncrementOld;
-				if (Theta < 0.99) {
+				if(ANY(Theta >= 0.99)) //! Non-convergence of Newton: Theta too large
+					break;
+				else
 					NewtonRate = Theta / (ONE - Theta);
-				}
-				else { //! Non-convergence of Newton: Theta too large
-					break;
-				}
-				if (NewtonIter < NewtonMaxit) {
-					//Predict error at the end of Newton process 
-					double NewtonPredictedErr = (NewtonIncrement * pow(Theta, (NewtonMaxit - NewtonIter - 1))) / (ONE - Theta);
-					if (NewtonPredictedErr >= NewtonTol) {
-						//Non-convergence of Newton: predicted error too large
-						double Qnewton = fmin(10.0, NewtonPredictedErr / NewtonTol);
-	                    Fac = 0.8 * pow(Qnewton, -ONE/((double)(NewtonMaxit-NewtonIter)));
-	                    break;
-					}
+				//Predict error at the end of Newton process 
+				double NewtonPredictedErr = (NewtonIncrement * pow(Theta, (NewtonMaxit - NewtonIter - 1))) / (ONE - Theta);
+				if(ANY (NewtonPredictedErr >= NewtonTol)) {
+					//Non-convergence of Newton: predicted error too large
+					double Qnewton = fmin(10.0, NewtonPredictedErr / NewtonTol);
+                    Fac = 0.8 * pow(Qnewton, -ONE/((double)(NewtonMaxit-NewtonIter)));
+                    break;
 				}
 			}
 

	int Nconsecutive = 0;
	int Nsteps = 0;
	while (t + Roundoff < t_end) {
		if (!Reject) {
			dydt (t, pr, y, F0);
		}
		if (!SkipLU) { 
			//need to update Jac/LU
			if (!SkipJac) {
				eval_jacob (t, pr, y, A);
			}
			RK_Decomp(H, E1, E2, A, ipiv1, ipiv2, &info);
			if (info != 0) {
				Nconsecutive += 1;
				if (Nconsecutive >= 5)
				{
					//todo implement return codes
					y[0] = logf(-1);
					return;
				}


            NewtonDone = (NewtonRate * NewtonIncrement <= NewtonTol);
#ifndef NEWTON_UNROLL
            if (ALL (NewtonDone)) break;
#else //only break if it's at the end of the unroll
            if(ALL (NewtonDone && (NewtonIter + 1) % NewtonIter == 0)) break;
#endif
            if (NewtonIter >= NewtonMaxit)
            {
				y[0] = logf(-1);
				return;
			}
		}
		if(ANY (!NewtonDone)) {
			H = Fac * H;
			Reject = true;
			SkipJac = true;
