diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 23263ae..76b2bb2 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -89,8 +89,9 @@ int uk(
     } else if (algorithm==ALGORITHM_RK4) {
       ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
     } else if (algorithm==ALGORITHM_RKF45) {
-      ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible);
+      ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, true);
     } else if (algorithm==ALGORITHM_RKBS23) {
+      // only compute k1 on the first step
       ns_step_rkbs23(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, time==starting_time);
     } else {
       fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
@@ -194,8 +195,9 @@ int enstrophy(
     } else if (algorithm==ALGORITHM_RK4) {
       ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
     } else if (algorithm==ALGORITHM_RKF45) {
-      ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible);
+      ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, true);
     } else if (algorithm==ALGORITHM_RKBS23) {
+      // only compute k1 on the first step
       ns_step_rkbs23(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, time==starting_time);
     } else {
       fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
@@ -348,8 +350,9 @@ int quiet(
     } else if (algorithm==ALGORITHM_RK4) {
       ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
     } else if (algorithm==ALGORITHM_RKF45) {
-      ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible);
+      ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, true);
     } else if (algorithm==ALGORITHM_RKBS23) {
+      // only compute k1 on the first step
       ns_step_rkbs23(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, time==starting_time);
     } else {
       fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
@@ -649,13 +652,17 @@ int ns_step_rkf45(
   _Complex double* k5,
   _Complex double* k6,
   _Complex double* tmp,
-  bool irreversible
+  bool irreversible,
+  // whether to compute k1 or leave it as is
+  bool compute_k1
 ){
   int kx,ky;
   double err,relative;
 
   // k1: u(t)
-  ns_rhs(k1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
+  if(compute_k1){
+    ns_rhs(k1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
+  }
 
   // k2 : u(t+1/4*delta)
   for(kx=0;kx<=K1;kx++){
@@ -724,7 +731,8 @@ int ns_step_rkf45(
   // error too big: repeat with smaller step
   else{
     *delta=factor*(*delta)*pow(relative*tolerance/err,0.2);
-    ns_step_rkf45(u,tolerance,factor,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible);
+    // do not recompute k1
+    ns_step_rkf45(u,tolerance,factor,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible,false);
   }
 
   return 0;
diff --git a/src/navier-stokes.h b/src/navier-stokes.h
index 308f108..3f1f6dc 100644
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@@ -56,7 +56,7 @@ int ns_step_rk4( _Complex double* u, int K1, int K2, int N1, int N2, double nu,
 // RK2
 int ns_step_rk2( _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2,fft_vect ifft, _Complex double* tmp1, _Complex double *tmp2, bool irreversible);
 // adaptive RK algorithm (Runge-Kutta-Fehlberg)
-int ns_step_rkf45( _Complex double* u, double tolerance, double factor, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double* k1, _Complex double* k2, _Complex double* k3, _Complex double* k4, _Complex double* k5, _Complex double* k6, _Complex double* tmp, bool irreversible);
+int ns_step_rkf45( _Complex double* u, double tolerance, double factor, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double* k1, _Complex double* k2, _Complex double* k3, _Complex double* k4, _Complex double* k5, _Complex double* k6, _Complex double* tmp, bool irreversible, bool compute_k1);
 // Runge-Kutta-Bogacki-Shampine
 int ns_step_rkbs23( _Complex double* u, double tolerance, double factor, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double** k1, _Complex double* k2, _Complex double* k3, _Complex double** k4, _Complex double* tmp, bool irreversible, bool compute_k1);