def dxmax_control(dxmax, dxmax0, j, j0, z, cycl, nr_success, b, n, nfev,
    minfac=1e-06, maxfac=10, reduce=8.0, increase=0.5, verbose=2):
    """Automatic control for incremental limit of system vector x. Based on the
    success of the nonlinear solution process the incremental limit will be
    either increased or decreased. If the newton iterations in the current inc.
    converged and the last 3 increments were successful (without any recycle)
    an increase of `dxmax` is performed. The amount of successful increments
    before the current one is read by the parameter `b`. The increase factor is
    based on a constant value `increase` and on the convergence rate of the
    current increment. The total increase factor is calculated as
    `increase_total = 1 + (nfev-n)/nfev * increase`. On the other hand, 
    if the current increment did not converge, `dxmax` will be reduced by a
    given value `decrease`. Both an increase and decrease is only performed
    inside the interval `minfac*dxmax0 <= dxmax <= maxfac*dxmax0`.
    
    Parameter
    ---------
    dxmax : ndarray
        current incremental limit of system vector x
    dxmax0 : ndarray
        initial incremental limit of system vector x
    j : int
        current control component of x
    j0 : int
        initial control component of x (at the beginning of the inc.)
    z : int
        number of recycles in current increment
    cycl: int
        max. number of recycles
    nr_success : bool
        flag to check if nonlinear solution process converged
    b : int
        number of converged increments before this increment
    n : int
        number of newton iterations until convergence was reached
    nfev : int
        allowed number of newton iterations
    minfac : float, optional
        minimum factor of dxmax0 where a decrease of dxmax is performed
        (`minfac*dxmax0 <= dxmax <= maxfac*dxmax0`)
    maxfac : float, optional
        minimum factor of dxmax0 where a decrease of dxmax is performed
        (`minfac*dxmax0 <= dxmax <= maxfac*dxmax0`)
    reduce : flaot, optional
        factor do decrease dxmax if current increment did not converge.
        (`dxmax= dxmax/reduce`)
    increase : float, optional
        factor do increase dxmax if current increment and the 3 increments
        before did converge.
        (`dxmax = dxmax * (1 + (nfev-n)/nfev * increase)`)
    verbose : int, optional
        Level of information during iterations (default is 0): 
        `verbose=0,1` ... no information, 
        `verbose=2`   ... print dxmax changes
    
    Returns
    -------
    dxmax : ndarray
        updated incremental limit of system vector x
    z : int
        updated number of recycles in current increment
    b : int
        updated number of converged increments before this increment
    dx_changed : bool
        flag to indicate if a change in dxmax was performed
    
    """
    dx_changed = False
    if not nr_success:
        if dxmax[0] > dxmax0[0] * minfac:
            if z + 1 == cycl:
                z = 0
                final_step_reduction = True
                if verbose > 1:
                    print(
                        '* reduce NR-step although control comp. changes (max. recycles reached).'
                        )
            else:
                final_step_reduction = False
                if verbose > 1:
                    print(
                        '* reduce NR-step size by factor: (inactive) due to recycle.'
                        )
            if j0 == j or final_step_reduction:
                if verbose > 1:
                    print('* reduce NR-step size by factor: {:10.3g}'.
                        format(1 / reduce))
                dxmax = dxmax / reduce
                dx_changed = True
                b = 0
    elif dxmax[0] < dxmax0[0] * maxfac:
        if j0 == j:
            if z > 0:
                pass
            elif b <= 3:
                if verbose > 1:
                    print(
                        '* increase NR-step (inactive). incs since last reduction: {0:1d}/{1:1d}'
                        .format(b, 3))
            else:
                dxmax = dxmax * (1 + (nfev - n) / nfev * increase)
                dx_changed = True
                if verbose > 1:
                    print('* increase NR-step size by factor: {:10.3g}'.
                        format(1 + (nfev - n) / nfev * increase))
    elif j == j0:
        dxmax = dxmax0.copy() * maxfac
        dx_changed = True
    return dxmax, z, b, dx_changed
