Numerical calculation of out-of-time order correlators (OTOCs)

Question

I want to numerically calculate OTOCs for different quantum mechanical systems. Consider the following Hamiltonian $$H=p_x^2+p_y^2+x^2y^2$$

and I want to calculate the following OTOC $$C_T(t)=-\left<[x(t),p_x(0)]^2\right>_\beta$$ where $$\left<\mathcal{O}\right>_\beta=Z^{-1}\text{Tr}(\mathcal{O} e^{-\beta H}),\qquad Z=\text{Tr}(e^{-\beta H})$$

You can look into the details for the recipe for calculating the OTOC in this paper: "Out-of-time-order correlators in quantum mechanics", in section 2.

My question is that I can only manage to calculate OTOC's using Mathematica, but it is very slow with numerical calculations. The essential idea behind the calculation of OTOCs is to find the energy eigenfunctions and energy eigenvalues for the given Hamiltonian. Then using the OTOC recipe, we can easily calculate the OTOC. Therefore I am looking for a programming language that can most efficiently calculate the eigenfunctions and energy eigenvalues of a given Hamiltonian, then calculate the matrix elements of the position operator $x_{nm}=\left<n\right|x\left|m\right>$ using parrallel computation. Parrallel computation is key since the matrix elements of the position operator form a quite large matrix if we want accurate results.

Do you know what programming language the following papers use

1. Out-of-time-order correlators in quantum mechanics
2. Lyapunov Exponent and Out-of-Time-Ordered Correlator's Growth Rate in a Chaotic System
3. Does scrambling equal chaos?
4. Gauging classical and quantum integrability through out-of-time ordered correlators