This is an amazing piece of code and I'm delighted to have found it! I think it might be useful for an application I am working on. The goal is to solve the chemical master equation

$$

\frac{\mathrm{d} P}{\mathrm{d} t} =\mathbf{A} P

$$

where @@\mathbf{A}@@ is a matrix of connections and @@P@@ represents the vector of probabilities of being in each state. (See this Wikipedia page for more details on master equations.) I got the idea for solving the master equation using the tensor formalism from this paper and its companion paper. Our system is similar in form to the CO oxidation model in the papers, but it has different chemical species and more complicated reactions.

To model this, I *think* I can use ITensor by applying an MPO (which may act on more than two sites!) to an MPS repeatedly and, at each timestep, measure the properties we are interested in (say, total number of atoms/molecules of each type). I haven't seen an example of this being done in ITensor yet, but this Wikipedia page suggests that a reaction-diffusion system could be modeled as a time-evolving set of non-quantum states.

**Has anyone done something like this in ITensor before? Are there any fundamental problems with using the library this way, particularly since this isn't exactly the canonical use case?**

I also understand that `toExpH`

doesn't work for MPOs that operate on more than two sites -- which, sadly, mine probably does. **Barring the obvious numerical instabilities, is there any reason why I couldn't do something like an explicit Euler integration instead of using the matrix exponential form to solve the above equation?**