def METHOD_NAME(mesh, lx, diff_flux, **model_options):
    r"""
    Consider a diffusion problem with a inhomogeneous Neumann condition and zero initial condition:
    .. math::
       c_t = \nu c_{xx}, c_x(0, t) = D, c_x(l, t) = 0, c(x, 0) = 0
    where :math:`D` is the diffusive flux boundary condition imposed, `diff_flux`.
    In order to solve it analytically, we decompose it into two diffusion problems:
     - a diffusion problem with homogeneous Neumann conditions and a nonzero initial condition,
    .. math::
           z_t = \nu z_{xx}, z_x(0, t) = 0, z_x(l, t) = 0, z(x, 0) = -I;
     - and a diffusion problem with homogeneous Neumann conditions and a nonzero source term,
    .. math::
           w_t = \nu w_{xx} + S, w_x(0, t) = 0, w_x(l, t) = 0, w(x, 0) = 0.
    Here :math:`I = I(x)` is set as :math:`I(x) = \alpha(x) D`, where
    :math:`\alpha(x) = -\frac{(l - x)^2}{2l}` and :math:`l` is the horizontal length of the domain.
    The source term :math:`S = S(x, t)` is given by :math:`S = I_t - \nu I_xx = -\nu\frac Dl`.
    Solving the inhomogeneous problem amounts to summing the solutions of the homogeneous problems
    and subtracting :math:`I`.
    The exact solution takes the form of a Fourier series, which we truncate appropriately.
    """
    x, y = SpatialCoordinate(mesh)
    nu = model_options['horizontal_diffusivity']
    time = model_options['simulation_end_time']
    # Initial condition and source term for two homogeneous Neumann problems
    P1 = get_functionspace(mesh, 'CG', 1)
    ic = Function(P1).interpolate(diff_flux*0.5*(lx - x)*(lx - x)/lx)
    source = Constant(-nu*diff_flux/lx, domain=mesh)
    # The solution uses truncated Fourier expansions, meaning we need the following...
    def phi(n):
        return cos(n*pi*x/lx)
    def ic_fourier_coeff(n):
        return assemble(2/lx*ic*phi(n)*dx)
    def source_fourier_coeff(n):
        return assemble(2/lx*source*phi(n)*dx)
    def source_term(n):
        """
        A simple quadrature routine is used to approximate the time integral.
        """
        I = 0
        tau = 0
        dt = 0.05
        while tau < time - 0.5*dt:
            I += exp(-nu*(n*pi/lx)**2*(t-tau))
            tau += dt
        I *= source_fourier_coeff(n)
        return I*phi(n)
    def ic_term(n):
        return ic_fourier_coeff(n)*exp(-nu*(n*pi/lx)**2*time)*phi(n)
    # Assemble truncated Fourier expansion
    sol = Function(P1, name='Fourier expansion')
    num_terms_source = 1  # Only one needed since source is constant
    num_terms_ic = 100
    expr = Constant(0.5*source_fourier_coeff(0)*time)
    expr = expr + Constant(0.5*ic_fourier_coeff(0))
    for k in range(1, num_terms_source):
        expr = expr + source_term(k)
    for k in range(1, num_terms_ic):
        expr = expr + ic_term(k)
    expr -= ic
    sol.interpolate(-expr)
    return sol