docs: 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