The below (< 50 lines) of code is the simplest way I could illustrate my possible misunderstanding. To be clear, this is *not* the actual problem I am trying to solve in my research (the solution is trivial); I'm just trying to find my way around ITensor and MPO/MPS.

Here is my current understanding of what this code **should** do: I construct an MPS of all "occupied" sites and two MPOs. The first is the Hamiltonian, which is supposed to subtract particles from occupied sites to yield empty sites, at a simple rate of 1. At each timestep, I apply the Hamiltonian times a small timestep @@h@@. This is the explicit Euler solution to the equation:

$$

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

$$

The observable I want to record is the number of occupied sites, so I use the second MPO, constructed from the number operator, and apply it as @@\langle \psi, B \psi \rangle@@.

Examining the raw output, or even plotting the results, is not encouraging, because the total number of occupied states is increasing, not decreasing. I'm sure I've just misunderstood something along the way, but I can't think of any particcular reason why the method I've described wouldn't yield (something close to) the analytic solution,

$$

y(t) = 3 e^{-t}.

$$

As always, help is greatly appreciated.

```
#include <itensor/all_basic.h>
#include <itensor/all_mps.h>
#include <iostream>
int main(int argc, char *argv[])
{
int N = 3;
auto sites = itensor::Spinless(N);
auto ampo = itensor::AutoMPO(sites);
auto ampo2 = itensor::AutoMPO(sites);
auto state = itensor::InitState(sites);
for (auto j : itensor::range1(N))
{
ampo += "A", j;
ampo2 += "N", j;
state.set(j, "Occ");
}
auto mpo = itensor::toMPO<itensor::ITensor>(ampo, { "Exact=", true });
auto mpo2 = itensor::toMPO<itensor::ITensor>(ampo2, { "Exact=", true });
auto mps = itensor::MPS(state);
double h = 1.e-6, t = 0;
for (auto i : itensor::range(1000))
{
mps.plusEq(h * itensor::applyMPO(mpo, mps, { "Exact=", true }));
t += h;
auto norm = itensor::overlap(mps, mps);
auto count = itensor::overlap(mps, mpo2, mps);
std::cout << std::scientific << t << " " << norm << " " << count
<< std::defaultfloat << std::endl;
}
return 0;
}
```