diff --git a/README.md b/README.md
index 8e866ae..97150f0 100644
--- a/README.md
+++ b/README.md
@@ -93,8 +93,9 @@ should be a `;` sperated list of `key=value` pairs. The possible keys are
 
 * `random_seed` (int, default ): seed for random initialization.
 
-* `algorithm`: either `RK4` for Runge-Kutta 4, `RK2` for Runge-Kutta 2, or
-  `RKF45` for the Runge-Kutta-Fehlberg adaptive step method.
+* `algorithm`: `RK4` for Runge-Kutta 4, `RK2` for Runge-Kutta 2, `RKF45` for
+  the Runge-Kutta-Fehlberg adaptive step method, `RKBS23` for the
+  Runge-Kutta-Bogacki-Shampine method.
 
 * `adaptive_tolerance` (double, default 1e-11): when using an adaptive step
   method, this is the maximal allowed relative error.
diff --git a/src/constants.cpp b/src/constants.cpp
index 8bd85ed..70a3e8b 100644
--- a/src/constants.cpp
+++ b/src/constants.cpp
@@ -34,4 +34,5 @@ limitations under the License.
 #define ALGORITHM_ADAPTIVE_THRESHOLD 1000
 // adaptive algorithms: index > ALGORITHM_ADAPTIVE_THRESHOLD
 #define ALGORITHM_RKF45 1001
+#define ALGORITHM_RKBS23 1002
 
diff --git a/src/main.c b/src/main.c
index 8cea261..dadba17 100644
--- a/src/main.c
+++ b/src/main.c
@@ -265,6 +265,9 @@ int print_params(
   case ALGORITHM_RKF45:
     fprintf(file,", algorithm=RKF45, tolerance=%.15e, factor=%.15e",parameters.adaptive_tolerance, parameters.adaptive_factor);
     break;
+  case ALGORITHM_RKBS23:
+    fprintf(file,", algorithm=RKBS23, tolerance=%.15e, factor=%.15e",parameters.adaptive_tolerance, parameters.adaptive_factor);
+    break;
   default:
     fprintf(file,", algorithm=RK4");
     break;
@@ -671,6 +674,9 @@ int set_parameter(
     else if (strcmp(rhs,"RKF45")==0){
       parameters->algorithm=ALGORITHM_RKF45;
     }
+    else if (strcmp(rhs,"RKBS23")==0){
+      parameters->algorithm=ALGORITHM_RKBS23;
+    }
     else{
       fprintf(stderr, "error: unrecognized algorithm '%s'\n",rhs);
       return(-1);
diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 811d4b8..23263ae 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -90,6 +90,8 @@ int uk(
       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);
+    } else if (algorithm==ALGORITHM_RKBS23) {
+      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);
     }
@@ -193,6 +195,8 @@ int enstrophy(
       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);
+    } else if (algorithm==ALGORITHM_RKBS23) {
+      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);
     }
@@ -216,7 +220,7 @@ int enstrophy(
       avg_en_x_a*=print_freq/(time-prevtime);
 
       // print to stderr so user can follow along
-      if(algorithm==ALGORITHM_RKF45){
+      if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
 	fprintf(stderr,"% .8e % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",time, avg_a, avg_en, avg_en_x_a, alpha, enstrophy, alpha*enstrophy, step);
 	printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy, step);
       } else {
@@ -345,6 +349,8 @@ int quiet(
       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);
+    } else if (algorithm==ALGORITHM_RKBS23) {
+      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);
     }
@@ -400,6 +406,12 @@ int ns_init_tmps(
     *tmp5=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
     *tmp6=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
     *tmp7=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
+  } else if (algorithm==ALGORITHM_RKBS23){
+    *tmp1=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
+    *tmp2=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
+    *tmp3=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
+    *tmp4=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
+    *tmp5=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
   } else {
     fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
   };
@@ -459,6 +471,12 @@ int ns_free_tmps(
     free(tmp5);
     free(tmp6);
     free(tmp7);
+  } else if (algorithm==ALGORITHM_RKBS23){
+    free(tmp1);
+    free(tmp2);
+    free(tmp3);
+    free(tmp4);
+    free(tmp5);
   } else {
     fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
   };
@@ -712,6 +730,106 @@ int ns_step_rkf45(
   return 0;
 }
 
+
+// next time step
+// adaptive RK algorithm (Runge-Kutta-Bogacki-Shampine method)
+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,
+  // the pointers k1 and k4 will be exchanged at the end of the routine
+  _Complex double** k1,
+  _Complex double* k2,
+  _Complex double* k3,
+  _Complex double** k4,
+  _Complex double* tmp,
+  bool irreversible,
+  // whether to compute k1
+  bool compute_k1
+){
+  int kx,ky;
+  double err,relative;
+
+  // k1: u(t)
+  // only compute it if it is the first step (otherwise, it has already been computed due to the First Same As Last property)
+  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++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)/2*(*k1)[klookup_sym(kx,ky,K2)];
+    }
+  }
+  ns_rhs(k2, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
+
+  // k3 : u(t+3/4*delta)
+  for(kx=0;kx<=K1;kx++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(3./4*k2[klookup_sym(kx,ky,K2)]);
+    }
+  }
+  ns_rhs(k3, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
+
+  // k4 : u(t+delta)
+  // tmp cpmputed here is the next step
+  for(kx=0;kx<=K1;kx++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(2./9*(*k1)[klookup_sym(kx,ky,K2)]+1./3*k2[klookup_sym(kx,ky,K2)]+4./9*k3[klookup_sym(kx,ky,K2)]);
+    }
+  }
+  ns_rhs(*k4, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
+  
+  // compute error
+  err=0;
+  relative=0;
+  for(kx=0;kx<=K1;kx++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+      // difference between 5th order and 4th order
+      err+=cabs((*delta)*(5./72*(*k1)[klookup_sym(kx,ky,K2)]-1./12*k2[klookup_sym(kx,ky,K2)]-1./9*k3[klookup_sym(kx,ky,K2)]+1./8*(*k4)[klookup_sym(kx,ky,K2)]));
+      relative+=cabs(tmp[klookup_sym(kx,ky,K2)]-u[klookup_sym(kx,ky,K2)]);
+    }
+  }
+
+  // compare relative error with tolerance
+  if(err<relative*tolerance){
+    // add to output
+    for(kx=0;kx<=K1;kx++){
+      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+	u[klookup_sym(kx,ky,K2)]=tmp[klookup_sym(kx,ky,K2)];
+      }
+    }
+    // next delta to use in future steps
+    *next_delta=(*delta)*pow(relative*tolerance/err,1./3);
+
+    // k1 in the next step will be this k4 (first same as last)
+    tmp=*k1;
+    *k1=*k4;
+    *k4=tmp;
+  }
+  // error too big: repeat with smaller step
+  else{
+    *delta=factor*(*delta)*pow(relative*tolerance/err,1./3);
+    // this will reuse the same k1 without re-computing it
+    ns_step_rkbs23(u,tolerance,factor,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,tmp,irreversible,false);
+  }
+
+  return 0;
+}
+
 // right side of Irreversible/Reversible Navier-Stokes equation
 int ns_rhs(
   _Complex double* out,
diff --git a/src/navier-stokes.h b/src/navier-stokes.h
index eff95fb..308f108 100644
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@@ -57,6 +57,8 @@ int ns_step_rk4( _Complex double* u, int K1, int K2, int N1, int N2, double nu,
 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);
+// 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);
 
 // right side of Irreversible/reversible Navier-Stokes equation
 int ns_rhs( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, bool irreversible);