 		dydt (t, pr, y, fy);
 
 		// Jacobian matrix
-		double A[NN * NN] = {ZERO};
 		eval_jacob (t, pr, y, A);
     	double gy[NN];
     	//gy = fy - A * y

#endif

	double beta = 0;
	while ((t < t_end) && (t + h > t)) {
		//initial krylov subspace sizes
		int m, m1, m2;

			
			if (err <= ONE) {

				#pragma unroll
				for (int i = 0; i < NN; ++i)
					sc[i] = f_temp[i];
				//memcpy(sc, f_temp, NN * sizeof(double));
				
				// minimum of classical and Gustafsson step size prediction
				h_new = fmin(h_new, (h / h_old) * pow((err_old / (err * err)), (1.0 / ORD)));
