diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 9799fad..ef03bf3 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -18,7 +18,6 @@ limitations under the License.
 #include "io.h"
 #include "navier-stokes.h"
 #include "complex_tools.h"
-#include <cblas.h>
 #include <math.h>
 #include <stdint.h>
 #include <stdlib.h>
@@ -90,21 +89,8 @@ int uk(
   // iterate
   time=starting_time;
   while(final_time<0 || time<final_time){
-    if(algorithm==ALGORITHM_RK2){
-      ns_step_rk2(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible, keep_en_cst, target_en);
-    } else if (algorithm==ALGORITHM_RK4) {
-      ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible, keep_en_cst, target_en);
-    } else if (algorithm==ALGORITHM_RKF45) {
-      ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, true);
-    } else if (algorithm==ALGORITHM_RKDP54) {
-      // only compute k1 on the first step
-      ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, time==starting_time);
-    } else if (algorithm==ALGORITHM_RKBS32) {
-      // only compute k1 on the first step
-      ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, keep_en_cst, target_en, time==starting_time);
-    } else {
-      fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
-    }
+    // step
+    ns_step(algorithm, u, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, starting_time, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en);
 
     time+=step;
     step=next_step;
@@ -206,21 +192,7 @@ int enstrophy(
   // iterate
   time=starting_time;
   while(final_time<0 || time<final_time){
-    if(algorithm==ALGORITHM_RK2){
-      ns_step_rk2(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible, keep_en_cst, target_en);
-    } else if (algorithm==ALGORITHM_RK4) {
-      ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible, keep_en_cst, target_en);
-    } else if (algorithm==ALGORITHM_RKF45) {
-      ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, true);
-    } else if (algorithm==ALGORITHM_RKDP54) {
-      // only compute k1 on the first step
-      ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, time==starting_time);
-    } else if (algorithm==ALGORITHM_RKBS32) {
-      // only compute k1 on the first step
-      ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, keep_en_cst, target_en, time==starting_time);
-    } else {
-      fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
-    }
+    ns_step(algorithm, u, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, starting_time, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en);
 
     time+=step;
 
@@ -387,21 +359,7 @@ int quiet(
   // iterate
   time=starting_time;
   while(final_time<0 || time<final_time){
-    if(algorithm==ALGORITHM_RK2){
-      ns_step_rk2(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible, keep_en_cst, target_en);
-    } else if (algorithm==ALGORITHM_RK4) {
-      ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible, keep_en_cst, target_en);
-    } else if (algorithm==ALGORITHM_RKF45) {
-      ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, true);
-    } else if (algorithm==ALGORITHM_RKDP54) {
-      // only compute k1 on the first step
-      ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, time==starting_time);
-    } else if (algorithm==ALGORITHM_RKBS32) {
-      // only compute k1 on the first step
-      ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, keep_en_cst, target_en, time==starting_time);
-    } else {
-      fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
-    }
+    ns_step(algorithm, u, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, starting_time, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en);
 
     time+=step;
     step=next_step;
@@ -554,6 +512,60 @@ int copy_u(
 }
 
 // next time step
+int ns_step(
+  unsigned int algorithm,
+  _Complex double* u,
+  int K1,
+  int K2,
+  int N1,
+  int N2,
+  double nu,
+  double* delta,
+  double* next_delta,
+  double adaptive_tolerance,
+  double adaptive_factor,
+  double max_delta,
+  unsigned int adaptive_norm,
+  double L,
+  _Complex double* g,
+  double time,
+  double starting_time,
+  fft_vect fft1,
+  fft_vect fft2,
+  fft_vect ifft,
+  // the pointers tmp1 and tmp2 will be exchanged at the end of the routine for first-same-as-last RK algorithms
+  _Complex double** tmp1,
+  _Complex double** tmp2,
+  _Complex double* tmp3,
+  _Complex double* tmp4,
+  _Complex double* tmp5,
+  _Complex double* tmp6,
+  _Complex double* tmp7,
+  bool irreversible,
+  bool keep_en_cst,
+  double target_en
+){
+  if(algorithm==ALGORITHM_RK2){
+    ns_step_rk2(u, K1, K2, N1, N2, nu, *delta, L, g, fft1, fft2, ifft, *tmp1, *tmp2, irreversible, keep_en_cst, target_en);
+  } else if (algorithm==ALGORITHM_RK4) {
+    ns_step_rk4(u, K1, K2, N1, N2, nu, *delta, L, g, fft1, fft2, ifft, *tmp1, *tmp2, tmp3, irreversible, keep_en_cst, target_en);
+  } else if (algorithm==ALGORITHM_RKF45) {
+    ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, delta, next_delta, L, g, fft1, fft2, ifft, *tmp1, *tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, true);
+  } else if (algorithm==ALGORITHM_RKDP54) {
+    // only compute k1 on the first step
+    // first-same-as-last with 2-nd argument
+    ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, delta, next_delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, time==starting_time);
+  } else if (algorithm==ALGORITHM_RKBS32) {
+    // only compute k1 on the first step
+    // first-same-as-last with 4-th argument
+    ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, delta, next_delta, L, g, fft1, fft2, ifft, tmp1, tmp3, tmp4, tmp2, tmp5, irreversible, keep_en_cst, target_en, time==starting_time);
+  } else {
+    fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers.edu\n",algorithm);
+  }
+
+  return (0);
+}
+
 // RK 4 algorithm
 int ns_step_rk4(
   _Complex double* u,
@@ -1343,7 +1355,7 @@ int ns_T_nofft(
   return 0;
 }
 
-
+/*
 // Jacobian of rhs
 int ns_D_rhs(
   _Complex double* out,
@@ -1365,7 +1377,6 @@ int ns_D_rhs(
     alpha=0;
   }
 
-  #define MATSIZE (K1*(2*(K2+1)+K2))
   for(i=0; i<MATSIZE*MATSIZE; i++){
     out[i]=0;
   }
@@ -1414,6 +1425,7 @@ _Complex double ns_d_T(
     return (ky*lx-kx*ly)*(qx*qx+qy*qy-lx*lx-ly*ly)/sqrt(lx*lx+ly*ly)/sqrt(qx*qx+qy*qy)*u[klookup_sym(-qx,-qy,K2)];
   }
 }
+*/
 
 
 // compute alpha
diff --git a/src/navier-stokes.h b/src/navier-stokes.h
index 37fd076..ae35bde 100644
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@@ -51,6 +51,7 @@ int ns_free_tmps( _Complex double* u, _Complex double* tmp1, _Complex double *tm
 int copy_u( _Complex double* u, _Complex double* u0, int K1, int K2);
 
 // next time step for Irreversible/reversible Navier-Stokes equation
+int ns_step( unsigned int algorithm, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, double L, _Complex double* g, double time, double starting_time, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double** tmp1, _Complex double** tmp2, _Complex double* tmp3, _Complex double* tmp4, _Complex double* tmp5, _Complex double* tmp6, _Complex double* tmp7, bool irreversible, bool keep_en_cst, double target_en);
 // RK4
 int ns_step_rk4( _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, _Complex double *tmp3, bool irreversible, bool keep_en_cst, double target_en);
 // RK2
@@ -71,10 +72,12 @@ int ns_T( _Complex double* u, int K1, int K2, int N1, int N2, fft_vect fft1, fft
 // convolution term in right side of equation (computed without fft)
 int ns_T_nofft( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2);
 
+/*
 // Jacobian of rhs
 int ns_D_rhs( _Complex double* out, _Complex double* u, int K1, int K2, double nu, double L, bool irreversible);
 // derivative of T with respect to u_k
 _Complex double ns_d_T( _Complex double* u, int kx, int ky, int lx, int ly, int K2);
+*/
 
 // compute alpha
 double compute_alpha( _Complex double* u, int K1, int K2, _Complex double* g, double L);