Floquet Quasi-Energy spectrum in periodically driven lattice (Numerics)

I am trying to understand the steps for numerically diagonalizing the Hamiltonian for an amplitude modulated lattice. I am trying to follow along with the calculations in section 4.6 of this thesis. The general idea is that we have our hamiltonian:
$$H = frac{p^2}{2m}+V_0cos^2(k_Lx)$$
I am familiar with solving this in the fourier basis, where $psi=e^{iqx}u_{nq}(x)$ and then $u_{nq}(x)$ is expanded into plane waves $u_{nq} = sum_m c_m e^{iQ_m x}$ to get the tridiagonal matrix which is diagonalized for a given q (example here in section 1.2.1.)

The modulated lattice becomes (making the d/dx substitution):
$$H = partial_x^2+V_0(1+alpha sin(omega t))cos^2(k_Lx)$$
From the reference (the above equation is 4.23, and is defined as h(t)):

Calculation of the quasienergy band structure is done by numerical
diagonalization of the Floquet Hamiltonian derived from the
stroboscopic time evolution operator U(T). This done by calculating
the time dependent Hamiltonian H({ti} , k) over a discrete mesh {ti}
with a time dependent amplitude of equation 4.23.

So the definitions are (choosing the exponential midpoint)
$$U(T)=prod_{n=1}^NU(t_n)$$
$$U(t+Delta t)=exp(-iDelta t h(t+Delta T/2))$$
$$G(k) propto log(U(T))$$
And the quasi energies are the eigenvalues of G(k).

My confusion comes into calculating the matrix U at a given time step. To calculate h(t), I need to pick a basis. Given that I calculate the time dependent bands in the plane wave representation, would it make sense to do that, where the above definition of $U(t_i)$ becomes the tridiagonal matrix as before, but with the lattice potential $V_0$ replaced by its time dependent value?

I tried coding this idea up in python to reproduce this (used the parameters from the paper with 400 time steps and 5 plane waves):

# Define Fourier Basis Hamiltonian
l = np.arange(-2,3)
def h(V0,q,l):
    diags = [V0 / 2 + Er*(q + 2 * k) ** 2 for k in l]
    Hmat = np.diag(diags)
    Hmat += V0 / 4 * (np.diag(np.ones(len(l) - 1), 1) + np.diag(np.ones(len(l) - 1), -1))
    return Hmat 

for q in qx:
    # Static
    Hmat = h(V0,q,l)
    evals, evecs = LA.eig(Hmat)
    Estatic.append(sorted(evals))
    # Initialize U
    Umatrix = np.eye(len(l))
    for t in times:
        # Get lattice potential at time = t
        v = V0*(1+am_amp*np.sin(am_freq*t))
        h1 = h(v,q,l)
        # U(t) = exp(-i delta t *h(tj) )
        Umatrix = np.matmul(LA.expm(-1j*dt*h1) , Umatrix)
    amEigs, amVecs = LA.eig(1j*np.log(Umatrix))

And my attempt at reproducing is nothing like figure 4.19 from the thesis. So I am curious if I have the right idea with the time dependent U calculation.