 	int Nconsecutive = 0;
 	int Nsteps = 0;
 	while (t + Roundoff < t_end) {
+		PRINT("Reject:%d\n", ANY(!Reject));
 		if(ANY(!Reject)) {
 			dydt (t, pr, y, F0);
 		}
+		PRINT("LU:%d\n", ANY(!SkipLU));
 		if(ANY(!SkipLU)) { 
 			//need to update Jac/LU
+			PRINT("Jac:%d\n", ANY(!SkipJac));
 			if(ANY(!SkipJac)) {
 				eval_jacob (t, pr, y, A);
 			}
 			RK_Decomp(H, E1, E2, A, ipiv1, ipiv2, &info);
+			PRINT("info:%d\n", ANY(!info));
 			if(ANY(info != 0)) {
 				Nconsecutive += 1;
 				if (Nconsecutive >= 5)

			if (NewtonIter > 0) 
			{
				Theta = NewtonIncrement / NewtonIncrementOld;
				if(ANY(Theta >= 0.99)) //! Non-convergence of Newton: Theta too large
					break;
				else
					NewtonRate = Theta / (ONE - Theta);
				//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);

				return;
			}
		}
		PRINT("NotDone:%d\n", ANY(!NewtonDone));
		if(ANY (!NewtonDone)) {
			H = Fac * H;
			Reject = true;
