def _make_diffusion_matrix(K, weight1=None, weight2=None):
    """Builds the general diffusion matrix with dimension nxn.

    .. note::

        :math:`n`   = number of points of diffusion axis
        :math:`n+1` = number of bounts of diffusion axis


    **Function-all argument** \n

    :param array K:         dimensionless diffusivities at cell boundaries
                            *(size: 1xn+1)*
    :param array weight1:   weight_1 *(size: 1xn+1)*
    :param array weight2:   weight_2 *(size: 1xn)*
    :returns:               completely listed tridiagonal diffusion matrix *(size: nxn)*
    :rtype:                 array

    .. note::

        The elements of array K are acutally dimensionless:

        .. math::

            K[i] = K_{\\textrm{physical}}  \\frac{\\Delta t}{(\\Delta y)^2}

        where :math:`K_{\\textrm{physical}}` is in unit :math:`\\frac{\\textrm{length}^2}{\\textrm{time}}`


    The diffusion matrix is build like the following

    .. math::

        \\textrm{diffTriDiag}=
        \\left[ \\begin{array}{cccccc}
        1+\\frac{s_1 }{w_{2,0}} & -\\frac{s_1}{w_{2,0}} & 0 &  & ... & 0  \\\\
        -\\frac{s_1}{w_{2,1}} & 1+\\frac{s_1 + s_2}{w_{2,1}} & -\\frac{s_2}{w_{2,1}} & 0 & ... & 0 \\\\
        0 & -\\frac{s_2}{w_{2,2}}  & 1+\\frac{s_2 + s_3}{w_{2,2}} & -\\frac{s_3}{w_{2,2}} &... & 0  \\\\
          &  & \\ddots & \\ddots & \\ddots & \\\\
        0 & 0 & ... & -\\frac{s_{n-2}}{w_{2,n-2}}  & 1+\\frac{s_{n-2} + s_{n-1}}{w_{2,{n-2}}} & -\\frac{s_{n-1}}{w_{2,{n-2}}} \\\\
        0 & 0 & ... & 0 & -\\frac{s_{n-1}}{w_{2,n-1}}  & 1+\\frac{s_{n-1}}{w_{2,n-1}} \\\\
        \\end{array} \\right]

    where

    .. math::

           \\begin{array}{lllllll}
                K   &= [K_0,     &K_1,    &K_2,    &...,&K_{n-1},  &K_{n}] \\\\
                w_1 &= [w_{1,0}, &w_{1,1},&w_{1,2},&...,&w_{1,n-1},&w_{1,n}] \\\\
                w_2 &= [w_{2,0}, &w_{2,1},&w_{2,2},&...,&w_{2,n-1}]
           \\end{array}

    and following subsitute:

    .. math::

        s_i = w_{1,i} K_i

    """

#           \\begin{eqnarray}
#              y    & = & ax^2 + bx + c \\\\
#              f(x) & = & x^2 + 2xy + y^2
#           \\end{eqnarray}

#    .. math::
#
#        K   &= [K_0,        &K_1,        &K_2,        &...    ,    &K_{n-1},      &K_{n}] \\\\
#        w_1 &= [w_{1,0},    &w_{1,1},    &w_{1,2},    &...    ,    &w_{1,n-1}, \\ &w_{1,n}] \\\\
#        w_2 &= [w_{2,0}, \\ &w_{2,1}, \\ &w_{2,2}, \\ &... \\ , \\ &w_{2,n-1}]    &o \\\\
#
#    """
    J = K.size - 1
    if weight1 is None:
        weight1 = np.ones_like(K)
    if weight2 is None:
        weight2 = np.ones(J)
    weightedK = weight1 * K
    Ka1 = weightedK[0:J] / weight2
    Ka3 = weightedK[1:J+1] / weight2
    Ka2 = np.insert(Ka1[1:J], 0, 0) + np.append(Ka3[0:J-1], 0)
    #  Atmosphere tridiagonal matrix
    #  this code makes a 3xN matrix, suitable for use with solve_banded
    #diag = np.empty((3, J))
    #diag[0, 1:] = -Ka3[0:J-1]
    #diag[1, :] = 1 + Ka2
    #diag[2, 0:J-1] = -Ka1[1:J]
    #  Build the full banded matrix instead
    A = (np.diag(1 + Ka2, k=0) +
         np.diag(-Ka3[0:J-1], k=1) +
         np.diag(-Ka1[1:J], k=-1))
    return A