CCTK_REAL state2 TYPE=gf TAGS='rhs="TestODESolvers::rhs2"' "State vector, stepped with twice the timestep"
CCTK_REAL rhs2 TYPE=gf "RHS of state vector"
CCTK_REAL error2 TYPE=gf "Error in state vector"
CCTK_REAL corder TYPE=gf "Estimate for convergence order"

STORAGE: state rhs

STORAGE: state rhs error state2 rhs2 error2 corder

  WRITES: state2(everywhere)

  READS: state
  WRITES: rhs

  READS: state, state2
  WRITES: rhs, rhs2

  SYNC: state

  SYNC: state, state2

  READS: state
  WRITES: error

  READS: state, state2
  WRITES: error, error2, corder

        state2[ind] = u0;

  // I want to step state2 with twice the stepsize, so need to adjust the time
  // at the RK substeps to have t_(k) = t0 + k dt * 2
  // TODO: this will break when using more than 1 thread
  static int last_it = -1;
  static CCTK_REAL original_time;
#pragma omp critical
  if(last_it != cctk_iteration) {
    last_it = cctk_iteration;
    original_time = cctk_time;
  }
  const CCTK_REAL cctk_time2 = original_time + 2*(cctk_time - original_time);

        rhs[ind] = order * pow(1 + cctk_time, order - 1);

        // solving u(t) for (1 + t)^(order - 1) = u_inverse gives:
        const CCTK_REAL u_inverse = pow(state[ind], (order - 1) / CCTK_REAL(order));
        // mix in some of the state vector and some of the analytic RHS
        rhs[ind] = 0.5*order*(u_inverse + pow(1. + cctk_time, order - 1));
        if(cctk_iteration % 2) {
          const CCTK_REAL u_inverse2 = pow(state2[ind], (order - 1) / CCTK_REAL(order));
          rhs2[ind] = 2*0.5*order*(u_inverse2 + pow(1. + cctk_time2, order - 1));
        } else {
          rhs2[ind] = 0.;
        }

        error2[ind] = state2[ind] - u0;
        if(error[ind] == 0 || error2[ind] == 0) { // happens reliably at cctk_time == cctk_initial_time
          corder[ind] = 0;
        } else {
          corder[ind] = log(abs(error2[ind]/error[ind]))/log(2.);
        }