def get_coeff_values(self):
    """
    Calculates the matrix of coefficients that is later used via sparse matrix
    solution techniques (scipy.sparse.linalg.spsolve) to compute the flexural
    response to the load. This step need only be performed once, and the
    coefficient matrix can very rapidly compute flexural solutions to any load.
    This makes this particularly good for probelms with time-variable loads or 
    that require iteration (e.g., water loading, in which additional water 
    causes subsidence, causes additional water detph, etc.).

    These must be linearly combined to solve the equation.

    13 coefficients: 13 matrices of the same size as the load
    
    NOTATION FOR COEFFICIENT BIULDING MATRICES (e.g., "cj0i_2"):
    c = "coefficient
    j = columns = x-value
    j0 = no offset: look at center of array
    i = rows = y-value
    i_2 = negative 2 offset (i2 = positive 2 offset)
    """

    # don't want to keep typing "self." everwhere!
    D = self.D
    drho = self.drho
    dx4 = self.dx4
    dy4 = self.dy4
    dx2dy2 = self.dx2dy2
    nu = self.nu
    g = self.g

    if np.isscalar(self.Te):
      # So much simpler with constant D! And symmetrical stencil
      self.cj2i0 = D/dy4
      self.cj1i_1 = 2*D/dx2dy2
      self.cj1i0 = -4*D/dy4 - 4*D/dx2dy2
      self.cj1i1 = 2*D/dx2dy2
      self.cj0i_2 = D/dx4
      self.cj0i_1 = -4*D/dx4 - 4*D/dx2dy2
      self.cj0i0 = 6*D/dx4 + 6*D/dy4 + 8*D/dx2dy2 + drho*g
      self.cj0i1 = -4*D/dx4 - 4*D/dx2dy2 # Symmetry
      self.cj0i2 = D/dx4 # Symmetry
      self.cj_1i_1 = 2*D/dx2dy2 # Symmetry
      self.cj_1i0 = -4*D/dy4 - 4*D/dx2dy2 # Symmetry
      self.cj_1i1 = 2*D/dx2dy2 # Symmetry
      self.cj_2i0 = D/dy4 # Symmetry
      # Bring up to size
      self.cj2i0 *= np.ones(self.qs.shape)
      self.cj1i_1 *= np.ones(self.qs.shape)
      self.cj1i0 *= np.ones(self.qs.shape)
      self.cj1i1 *= np.ones(self.qs.shape)
      self.cj0i_2 *= np.ones(self.qs.shape)
      self.cj0i_1 *= np.ones(self.qs.shape)
      self.cj0i0 *= np.ones(self.qs.shape)
      self.cj0i1 *= np.ones(self.qs.shape)
      self.cj0i2 *= np.ones(self.qs.shape)
      self.cj_1i_1 *= np.ones(self.qs.shape)
      self.cj_1i0 *= np.ones(self.qs.shape)
      self.cj_1i1 *= np.ones(self.qs.shape)
      self.cj_2i0 *= np.ones(self.qs.shape)
      # Create coefficient arrays to manage boundary conditions
      self.cj2i0_coeff_ij = self.cj2i0.copy()
      self.cj1i_1_coeff_ij = self.cj1i_1.copy()
      self.cj1i0_coeff_ij = self.cj1i0.copy()
      self.cj1i1_coeff_ij = self.cj1i1.copy()
      self.cj0i_2_coeff_ij = self.cj0i_2.copy()
      self.cj0i_1_coeff_ij = self.cj0i_1.copy()
      self.cj0i0_coeff_ij = self.cj0i0.copy()
      self.cj0i1_coeff_ij = self.cj0i1.copy()
      self.cj0i2_coeff_ij = self.cj0i2.copy()
      self.cj_1i_1_coeff_ij = self.cj_1i_1.copy()
      self.cj_1i0_coeff_ij = self.cj_1i0.copy()
      self.cj_1i1_coeff_ij = self.cj_1i1.copy()
      self.cj_2i0_coeff_ij = self.cj_2i0.copy()
      
    elif type(self.Te) == np.ndarray:
    
      #######################################################
      # GENERATE COEFFICIENT VALUES FOR EACH SOLUTION TYPE. #
      #    "vWC1994" IS THE BEST: LOOSEST ASSUMPTIONS.      #
      #        OTHERS HERE LARGELY FOR COMPARISON           #
      #######################################################
      
      # All derivatives here, to make reading the equations below easier
      D00 = D[1:-1,1:-1]
      D10 = D[1:-1,2:]
      D_10 = D[1:-1,:-2]
      D01 = D[2:,1:-1]
      D0_1 = D[:-2,1:-1]
      D11 = D[2:,2:]
      D_11 = D[2:,:-2]
      D1_1 = D[:-2,2:]
      D_1_1 = D[:-2,:-2]
      # Derivatives of D -- not including /(dx^a dy^b)
      D0  = D00
      Dx  = (-D_10 + D10)/2.
      Dy  = (-D0_1 + D01)/2.
      Dxx = (D_10 - 2.*D00 + D10)
      Dyy = (D0_1 - 2.*D00 + D01)
      Dxy = (D_1_1 - D_11 - D1_1 + D11)/4.
      
      if self.PlateSolutionType == 'vWC1994':
        # van Wees and Cloetingh (1994) solution, re-discretized by me
        # using a central difference approx. to 2nd order precision
        # NEW STENCIL
        # x = -2, y = 0
        self.cj_2i0_coeff_ij = (D0 - Dx) / dx4
        # x = 0, y = -2
        self.cj0i_2_coeff_ij = (D0 - Dy) / dy4
        # x = 0, y = 2
        self.cj0i2_coeff_ij = (D0 + Dy) / dy4
        # x = 2, y = 0
        self.cj2i0_coeff_ij = (D0 + Dx) / dx4
        # x = -1, y = -1
        self.cj_1i_1_coeff_ij = (2.*D0 - Dx - Dy + Dxy*(1-nu)/2.) / dx2dy2
        # x = -1, y = 1
        self.cj_1i1_coeff_ij = (2.*D0 - Dx + Dy - Dxy*(1-nu)/2.) / dx2dy2
        # x = 1, y = -1
        self.cj1i_1_coeff_ij = (2.*D0 + Dx - Dy - Dxy*(1-nu)/2.) / dx2dy2
        # x = 1, y = 1
        self.cj1i1_coeff_ij = (2.*D0 + Dx + Dy + Dxy*(1-nu)/2.) / dx2dy2
        # x = -1, y = 0
        self.cj_1i0_coeff_ij = (-4.*D0 + 2.*Dx + Dxx)/dx4 + (-4.*D0 + 2.*Dx + nu*Dyy)/dx2dy2
        # x = 0, y = -1
        self.cj0i_1_coeff_ij = (-4.*D0 + 2.*Dy + Dyy)/dy4 + (-4.*D0 + 2.*Dy + nu*Dxx)/dx2dy2
        # x = 0, y = 1
        self.cj0i1_coeff_ij = (-4.*D0 - 2.*Dy + Dyy)/dy4 + (-4.*D0 - 2.*Dy + nu*Dxx)/dx2dy2
        # x = 1, y = 0
        self.cj1i0_coeff_ij = (-4.*D0 - 2.*Dx + Dxx)/dx4 + (-4.*D0 - 2.*Dx + nu*Dyy)/dx2dy2
        # x = 0, y = 0
        self.cj0i0_coeff_ij = (6.*D0 - 2.*Dxx)/dx4 \
                     + (6.*D0 - 2.*Dyy)/dy4 \
                     + (8.*D0 - 2.*nu*Dxx - 2.*nu*Dyy)/dx2dy2 \
                     + drho*g
                     
      elif self.PlateSolutionType == 'G2009':
        # STENCIL FROM GOVERS ET AL. 2009 -- first-order differences
        # x is j and y is i b/c matrix row/column notation
        # Note that this breaks down with b.c.'s that place too much control 
        # on the solution -- harmonic wavetrains
        # x = -2, y = 0
        self.cj_2i0_coeff_ij = D_10/dx4
        # x = -1, y = -1
        self.cj_1i_1_coeff_ij = (D_10 + D0_1)/dx2dy2
        # x = -1, y = 0
        self.cj_1i0_coeff_ij = -2. * ( (D0_1 + D00)/dx2dy2 + (D00 + D_10)/dx4 )
        # x = -1, y = 1
        self.cj_1i1_coeff_ij = (D_10 + D01)/dx2dy2
        # x = 0, y = -2
        self.cj0i_2_coeff_ij = D0_1/dy4
        # x = 0, y = -1
        self.cj0i_1_coeff_ij = -2. * ( (D0_1 + D00)/dx2dy2 + (D00 + D0_1)/dy4)
        # x = 0, y = 0
        self.cj0i0_coeff_ij = (D10 + 4.*D00 + D_10)/dx4 + (D01 + 4.*D00 + D0_1)/dy4 + (8.*D00/dx2dy2) + drho*g
        # x = 0, y = 1
        self.cj0i1_coeff_ij = -2. * ( (D01 + D00)/dy4 + (D00 + D01)/dx2dy2 )
        # x = 0, y = 2
        self.cj0i2_coeff_ij = D0_1/dy4
        # x = 1, y = -1
        self.cj1i_1_coeff_ij = (D10+D0_1)/dx2dy2
        # x = 1, y = 0
        self.cj1i0_coeff_ij = -2. * ( (D10 + D00)/dx4 + (D10 + D00)/dx2dy2 )
        # x = 1, y = 1
        self.cj1i1_coeff_ij = (D10 + D01)/dx2dy2
        # x = 2, y = 0
        self.cj2i0_coeff_ij = D10/dx4
      else:
        sys.exit("Not an acceptable plate solution type. Please choose from:\n"+
                  "* vWC1994\n"+
                  "* G2009\n"+
                  "")
                  
      ################################################################
      # CREATE COEFFICIENT ARRAYS: PLAIN, WITH NO B.C.'S YET APPLIED #
      ################################################################
      # x = -2, y = 0
      self.cj_2i0 = self.cj_2i0_coeff_ij.copy()
      # x = -1, y = -1
      self.cj_1i_1 = self.cj_1i_1_coeff_ij.copy()
      # x = -1, y = 0
      self.cj_1i0 = self.cj_1i0_coeff_ij.copy()
      # x = -1, y = 1
      self.cj_1i1 = self.cj_1i1_coeff_ij.copy()
      # x = 0, y = -2
      self.cj0i_2 = self.cj0i_2_coeff_ij.copy()
      # x = 0, y = -1
      self.cj0i_1 = self.cj0i_1_coeff_ij.copy()
      # x = 0, y = 0
      self.cj0i0 = self.cj0i0_coeff_ij.copy()
      # x = 0, y = 1
      self.cj0i1 = self.cj0i1_coeff_ij.copy()
      # x = 0, y = 2
      self.cj0i2 = self.cj0i2_coeff_ij.copy()
      # x = 1, y = -1
      self.cj1i_1 = self.cj1i_1_coeff_ij.copy()
      # x = 1, y = 0
      self.cj1i0 = self.cj1i0_coeff_ij.copy()
      # x = 1, y = 1
      self.cj1i1 = self.cj1i1_coeff_ij.copy()
      # x = 2, y = 0
      self.cj2i0 = self.cj2i0_coeff_ij.copy()

    # Provide rows and columns in the 2D input to later functions
    self.ncolsx = self.cj0i0.shape[1]
    self.nrowsy = self.cj0i0.shape[0]
Calculates the matrix of coefficients that is later used via sparse matrix
    solution techniques (scipy.sparse.linalg.spsolve) to compute the flexural
    response to the load. This step need only be performed once, and the
    coefficient matrix can very rapidly compute flexural solutions to any load.
    This makes this particularly good for probelms with time-variable loads or 
    that require iteration (e.g., water loading, in which additional water 
    causes subsidence, causes additional water detph, etc.).

    These must be linearly combined to solve the equation.

    13 coefficients: 13 matrices of the same size as the load
    
    NOTATION FOR COEFFICIENT BIULDING MATRICES (e.g., "cj0i_2"):
    c = "coefficient
    j = columns = x-value
    j0 = no offset: look at center of array
    i = rows = y-value
    i_2 = negative 2 offset (i2 = positive 2 offset)