def fast_inhomo_sweep_time_evolution(Ep, epsilonp, gamma,
                                     omega_level, rm, xi, theta,
                                     inhomogeneity,
                                     semi_analytic=True,
                                     file_name=None, return_code=False):
    r"""Return a spectrum of time evolutions of the density matrix.

    We test a basic two-level system.

    >>> from fast.bloch import phase_transformation
    >>> from fast import PlaneWave, electric_field_amplitude_top, Atom
    >>> Ne = 2
    >>> Nl = 1
    >>> a0 = physical_constants["Bohr radius"][0]
    >>> rm = [np.array([[0, 0], [a0, 0]]),
    ...       np.array([[0, 0], [0, 0]]),
    ...       np.array([[0, 0], [0, 0]])]
    >>> xi = np.array([[[0, 1], [1, 0]]])
    >>> omega_level = [0, 2.4e15]
    >>> gamma21 = 2*np.pi*6e6
    >>> gamma = np.array([[0, -gamma21], [gamma21, 0]])
    >>> theta = phase_transformation(Ne, Nl, rm, xi)

    >>> Ep, omega_laser = define_laser_variables(Nl)
    >>> laser = PlaneWave(0, 0, 0, 0)
    >>> epsilonp = [laser.epsilonp]
    >>> k = [laser.k]
    >>> detuning_knob = [symbols("delta1", real=True)]

    A map to unfold the density matrix.
    >>> unfolding = Unfolding(Ne, True, True, True)

    >>> Eps = electric_field_amplitude_top(1e-3, 1e-3, 1, "SI")
    >>> Eps = [Eps]

    >>> t = np.linspace(0, 1e-6, 11)
    >>> rho0 = np.array([[1, 0], [0, 0]])
    >>> rho0 = unfolding(rho0)

    We define the Doppler broadening.
    >>> Nvz = 15
    >>> shape = [Nvz]
    >>> stds = [[-4, 4]]
    >>> T = 273.15+20
    >>> mass = Atom("Rb", 87).mass
    >>> aux = (shape, stds, T, mass, detuning_knob, k,
    ...    omega_level, xi, theta, unfolding, ["z", "x", "y"],
    ...    True)

    >>> doppler_effect = DopplerBroadening(*aux)

    We get a function for the frequency sweep of time evolution.
    >>> aux = (Ep, epsilonp, gamma,
    ...        omega_level, rm, xi, theta,
    ...        doppler_effect,
    ...        True,
    ...        "eqs")

    >>> inhomo_time_evolution = fast_inhomo_sweep_time_evolution(*aux)

    >>> amp = 1000e6*2*np.pi
    >>> Ndelta = 101
    >>> detuning_knobs = [[-amp, amp, Ndelta]]

    >>> deltas, rhot = inhomo_time_evolution(t, rho0, Eps, detuning_knobs)
    >>> print rhot.shape
    (101, 11, 15, 3)

    """
    # We unpack variables.
    if True:
        Nl = xi.shape[0]
    # We determine which arguments are constants.
    if True:
        try:
            Ep = np.array([complex(Ep[l]) for l in range(Nl)])
            variable_Ep = False
        except:
            variable_Ep = True

        try:
            epsilonp = [np.array([complex(epsilonp[l][i]) for i in range(3)])
                        for l in range(Nl)]
            variable_epsilonp = False
        except:
            variable_epsilonp = True
    # We obtain code for the time evolution.
    if True:
        detuning_knob = symbols("delta1:"+str(Nl))
        args = (Ep, epsilonp, detuning_knob, gamma, omega_level, rm, xi, theta,
                file_name, True)

        args = (Ep, epsilonp, detuning_knob, gamma, omega_level, rm, xi,
                theta, inhomogeneity, True, file_name, True)
        inhomo_time_evolution = fast_inhomo_time_evolution(*args)
        code = inhomo_time_evolution+"\n\n"
    # We establish the arguments of the output function.
    if True:
        code += "def inhomo_sweep_time_evolution(t, rho0, "
        if variable_Ep: code += "Ep, "
        if variable_epsilonp: code += "epsilonp, "
        code += "detuning_knob, "
        code += "inhomo_time_evolution=inhomo_time_evolution):\n"
        code += '    r"""A fast frequency sweep of the steady state."""\n'
    # Code to determine the sweep range.
    if True:
        code += """    sweepN = -1\n"""
        code += """    for i, delta in enumerate(detuning_knob):\n"""
        code += """        if hasattr(delta, "__getitem__"):\n"""
        code += """            sweepN = i\n"""
        code += """            delta0 = delta[0]\n"""
        code += """            deltaf = delta[1]\n"""
        code += """            Ndelta = delta[2]\n"""
        code += """            break\n\n"""
        code += """    if sweepN == -1:\n"""
        code += """        s = 'One of the detuning knobs '\n"""
        code += """        s += 'must be of the form '\n"""
        code += """        s += '(start, stop, Nsteps)'\n"""
        code += """        raise ValueError(s)\n\n"""
        code += """    deltas = np.linspace(delta0, deltaf, Ndelta)\n\n"""
    # We call time_evolution.
    if True:
        code += "    args = [[t, rho0, "
        if variable_Ep: code += "Ep, "
        if variable_epsilonp: code += "epsilonp, "
        code += """list(detuning_knob[:sweepN]) +\n"""
        code += """            [deltas[i]] +\n"""
        code += """            list(detuning_knob[sweepN+1:])]\n"""
        code += """          for i in range(Ndelta)]\n\n"""
        code += "    rho = np.array([inhomo_time_evolution(*argsi)\n"
        code += "                   for argsi in args])\n\n"
    # We finish the code.
    if True:
        code += "    return deltas, rho\n"
    # We write the code to file if provided, and execute it.
    if True:
        if file_name is not None:
            f = file(file_name+".py", "w")
            f.write(code)
            f.close()

        inhomo_sweep_time_evolution = code
        if not return_code:
            exec inhomo_sweep_time_evolution
    return inhomo_sweep_time_evolution