def nonlinear_system(model, initial_dr=None, maxit=10, tol=1e-8,  grid={}, distribution={}, verbose=True):

    '''
    Finds a global solution for ``model`` by solving one large system of equations
    using a simple newton algorithm.

    Parameters
    ----------
    model: NumericModel
        "dtcscc" model to be solved
    verbose: boolean
        if True, display iterations
    initial_dr: decision rule
        initial guess for the decision rule
    maxit:  int
        maximum number of iterationsd
    tol: tolerance criterium for successive approximations
    grid: grid options
    distribution: distribution options

    Returns
    -------
    decision rule :
        approximated solution
    '''


    if verbose:
        headline = '|{0:^4} | {1:10} | {2:8} |'
        headline = headline.format('N', ' Error', 'Time')
        stars = '-'*len(headline)
        print(stars)
        print(headline)
        print(stars)
        # format string for within loop
        fmt_str = '|{0:4} | {1:10.3e} | {2:8.3f} |'

    f = model.functions['arbitrage']
    g = model.functions['transition']
    p = model.calibration['parameters']

    distrib = model.get_distribution(**distribution)
    nodes, weights = distrib.discretize()

    approx = model.get_grid(**grid)
    ms = create_interpolator(approx, approx.interpolation)

    grid = ms.grid

    if initial_dr is None:
        dr = approximate_controls(model)
    else:
        dr = initial_dr

    ms.set_values(dr(grid))

    x = dr(grid)
    x0 = x.copy()

    it = 0
    err = 10


    a0 = x0.copy().reshape((x0.shape[0]*x0.shape[1],))
    a = a0.copy()

    while err > tol and it < maxit:

        it += 1
        t1 = time.time()

        r, da = residuals(f, g, grid, a.reshape(x0.shape), ms, nodes, weights, p, diff=True)[:2]
        r = r.flatten()

        err = abs(r).max()

        t2 = time.time()

        if verbose:
            print(fmt_str.format(it, err, t2-t1))

        if err > tol:
            a -= scipy.sparse.linalg.spsolve(da, r)

    if verbose:
        print(stars)

    return ms
Finds a global solution for ``model`` by solving one large system of equations
    using a simple newton algorithm.

    Parameters
    ----------
    model: NumericModel
        "dtcscc" model to be solved
    verbose: boolean
        if True, display iterations
    initial_dr: decision rule
        initial guess for the decision rule
    maxit:  int
        maximum number of iterationsd
    tol: tolerance criterium for successive approximations
    grid: grid options
    distribution: distribution options

    Returns
    -------
    decision rule :
        approximated solution