def Heun_ndt(f_S, Si, dt, T, **kwargs):
    """ Solve dS/dt = f_S(S, t), S(t=t_i)=Si, at t=T with n steps (n=T/dt), using the explicit Heun's method

    Parameters
    ----------
    f_S : function
        State function for given model sub-domain (canopy, unsaturated zone, saturated zone)
    Si : float
        State at time t=t_i
    dt : float or int
        Time step [days]
    T : float or int
        Time period [days]
    kwargs : dict
        *kwargs* are used to specify the additional parameters used by the the state function (f_S)

    Returns
    -------
    dS : float
        Integrated state at t=ndt with n=T/dt
    """

    n = int(T / dt)
    dS = Si
    for _ in range(n):
        K1 = f_S(dS, **kwargs)
        K2 = f_S((K1 * dt) + dS, **kwargs)
        dS = dS + (0.5 * dt * (K1 + K2))
    return dS