Numerical solution of driven infinite well with Floquet matrix I am trying to understand and replicate figure 3 of this paper.
The idea is that we have our box with a periodic drive
$$H(t) = \frac{p^2}{2m}-F_0x\cos(\omega t)$$
The author states

In order to calculate the truncated time-evolution matrix $U(T, 0)$ one takes each of the unperturbed box eigenstates as initial condition, $ψ_n(x, 0) = ϕ_n(x)$, and computes the states $ψ_n(x, T )$ resulting after one period, collecting these state vectors as columns of the monodromy matrix

I am having difficulty figuring this out.  I know the unperturbed energies and wavefunctions $|\psi(0)\rangle$.  I also know that $|\psi(t)\rangle = U(t,0)|\psi(0)\rangle$
with ($\hbar=1$)
$$U(t,0)=\exp\left[-i\int_0^tH \ dt\right]$$
This is where I am unsure of how to numerically calculate.   Usually, I would do just $\psi(t)=\psi(0)\exp[-iE_n t]$, but because of the periodicity, we actually are able to expand
$$U(t,0)=\sum_n |n\rangle\langle n| \exp[-i\epsilon_n T]$$
In quasi-energies, which I want to solve for, so I am not sure of the time dependance of H to integrate to solve for $|\psi(x,t)\rangle$.  Any suggestions on how to approach this are appreciated.
 A: Here there's no fancy math that will save you - you just need to explicitly solve the TDSE numerically for a bunch of different initial conditions.
The algorithm looks more or less like this:


*

*For each $m=1,\ldots,N_\mathrm{max}$:


*

*Formulate the Schrödinger equation for the vector $(\psi_n(t))_{n=1}^{N_\mathrm{max}}$ as a set of coupled ordinary differential equations,
$$i\hbar \partial_t \psi_n(t) = E_n\psi_n(t)+\sum_{k=1}^{N_\mathrm{max}} Fx_{nk}\psi_k(t) \quad\text{under}\quad \psi_n(t) = \delta_{mn}$$

*Solve this set of coupled ODEs numerically up to $T=2\pi/\omega$.

*Store the results at the end of the period as $U_{mn} = \psi_n(T)$. This is known as the monodromy matrix. 


*Assemble the results of all the different $m$ into one big matrix $U_{mn}$.

*Diagonalize this matrix. Its eigenvalues are of the form $e^{i \varepsilon_k T/\hbar}$, i.e. the exponentiated quasienergies.

*Repeat for a different coupling.


This is what the authors mean when they say

then [we] have to deal with a system of $n_\mathrm{max}$ complex coupled ordinary differential equations which can be integrated numerically by standard routines.

I gave this a brief run in Mathematica and it's not that hard (though of course this gets harder and harder as $n_\mathrm{max}$ gets bigger). The results look roughly like this:

I haven't checked the constants for normalization or the convergence with respect to $n_\mathrm{max}$, so take that as a rough guide only; the code is at

Import["http://halirutan.github.io/Mathematica-SE-Tools/decode.m"]["http://i.stack.imgur.com/rPjH2.png"]

