For steady-state solutions, I have found that the most computationally efficient method is to reinterpret the kinetics network as a Markov process, calculate the transition-rate matrix, and then factorize it (using a special positivity-preserving factorization: https://arxiv.org/abs/2101.11657).
I have a PR for jaff that creates this matrix and does C++ codegen here that I am happy to port to jaco if this is of interest: jaff-chemistry/jaff#38
For steady-state solutions, I have found that the most computationally efficient method is to reinterpret the kinetics network as a Markov process, calculate the transition-rate matrix, and then factorize it (using a special positivity-preserving factorization: https://arxiv.org/abs/2101.11657).
I have a PR for jaff that creates this matrix and does C++ codegen here that I am happy to port to jaco if this is of interest: jaff-chemistry/jaff#38