I am trying to simulate a quantum system with the Bose-Hubbard model.

I have built my Hamiltonian, the initial state, the number operator matrices and everything is expressed in the same basis, so in principle there should be no problem here. Yet, when trying to time-evolve my initial state, things start to go wrong. For example, a lattice with size 5 (space for 5 bosons) and 2 bosons in positions 01000 and 00010 should give a simmetric evolution over time, since there is no side preference. But it’s not working for me, one side gets more particles and so on.

I have been looking at this problem for a few days now, but solved nothing thus far. As far as I’m concerned, the theoretical-physical aspect of my program is fine. There may, be some algebraic problem going on when multiplying the matrices and the vectors, but everything I have tried failed to solve the problem.

I leave my program here, please feel free to correct any kind of problem you see. Thanks a lot before hand.

https://github.com/mikelpuy/Bose-Hubbard/blob/main/2%20particles%20new.py