def METHOD_NAME(system, ii, jj, f):
    # Can be implemented as 2 calls of ``L``, but would use 2 ``dt`` which
    # cannot be different.
    def _func(t, rho, dt=1e-6):
        f0 = f(t, rho)
        bi = system.diffusion(t, rho)[ii]
        bj = system.diffusion(t, rho)[jj]
        out = rho *0.
        n = rho.shape[0]
        for i, j in product(range(n), repeat=2):
            dxi = basis(n, i, dtype="Dense").data
            dxj = basis(n, j, dtype="Dense").data
            sec = f(t, (rho + dxi * dt + dxj * dt))
            sec = sec - f(t, (rho + dxj * dt))
            sec = sec - f(t, (rho + dxi * dt))
            sec = sec + f0
            sec = sec / dt / dt
            out = out + (
                sec * _data.inner(dxi, bi) * _data.inner(dxj, bj)
            )
            df = (f(t, (rho + dxj * dt)) - f0) / dt
            db = (
                system.diffusion(t, (rho + dxi * dt))[jj]
                - system.diffusion(t, rho)[jj]
            ) / dt
            out = out + (
                df * _data.inner(dxi, bi) * _data.inner(dxj, db)
            )
        return out
    return _func