The problem wasn't actually with the RK method, but the erroneous test for probability values for a specific aa state. The appropriate test was whether the sum of probability values across aa state was 0. I have not checked any other ode solvers, but ode45 gives a 20 fold speed up over lsoda. As a result, I recommend checking the other solvers and updating ode routines to use the fastest one.
The problem wasn't actually with the RK method, but the erroneous test for probability values for a specific aa state. The appropriate test was whether the sum of probability values across aa state was 0. I have not checked any other ode solvers, but ode45 gives a 20 fold speed up over lsoda. As a result, I recommend checking the other solvers and updating ode routines to use the fastest one.