I’m trying to calculte matrix elements of teh Hamiltonian for teh harmonic oscillator. I use the solutions from sympy.physics.quatnum.qho_1d. My expectations would be that
<pis(i)|H|psi(j)>=0
if $i!=j$ where psi(i) is teh i’th solution of the harominc oscillator.
However, in my code, this is not the case if i+2=j. For all other combinations, it works.

import sympy as sp
from sympy.physics.qho_1d import psi_n
import sympy.physics.units
import numpy as np

x =sp.Symbol('x',real=True)
m=1
omega=1
hbar=sp.physics.units.hbar

def H(psi):
    return -hbar**2/(2*m) *sp.diff(psi,x,x) +1/2*m**2*omega**2*x**2*psi

psi = sp.Function('psi')


# computes matrix elements <psi1|H|psi2>
def matrix_elem(H,psi1,psi2):
   return sp.simplify(sp.integrate(sp.simplify(psi1.conjugate()*H(psi2)),(x,-sp.oo,sp.oo)))


max_n=5
H_matrix = sp.Matrix(
  [[matrix_elem(H,psi_n(i,x,1,1),psi_n(j,x,1,1)) for j in range(max_n)] for i in range(max_n)])

print(sp.simplify(H_matrix))
for i in range(max_n):
  for j in range(max_n):
    if i!=j and H_matrix[i,j] !=0: print(f"{i} {j}: {H_matrix[i,j]=}") 

An other thing that is not clear to me is that hbar**2/hbar is not simplified to hbar.

Many thanks in advance for any suggestions.