code: private static double[][] bcucof(double[] y, double[] y1, double[] y2, double[] y12, double d1, double d2) {
  
          double d1d2 = d1 * d2;
          double[] cl = new double[16];
          double[] x = new double[16];
          double[][] c = new double[4][4];
  
          for (int i = 0; i < 4; i++) {
              x[i] = y[i];
              x[i + 4] = y1[i] * d1;
              x[i + 8] = y2[i] * d2;
              x[i + 12] = y12[i] * d1d2;
          }
  
          double xx;
          for (int i = 0; i < 16; i++) {
              xx = 0.0;
              for (int k = 0; k < 16; k++) {
                  xx += wt[i][k] * x[k];
              }
              cl[i] = xx;
          }
  
          int l = 0;
          for (int i = 0; i < 4; i++) {
              for (int j = 0; j < 4; j++) {
                  c[i][j] = cl[l++];
              }
          }
  
          return c;
      }
docs: Given arrays y[0..3], y1[0..3], y2[0..3], and y12[0..3], containing
  the function, gradients, and cross-derivative at the four grid points
  of a rectangular grid cell (numbered counterclockwise from the lower
  left), and given d1 and d2, the length of the grid cell in the 1 and 2
  directions, returns the table c[0..3][0..3] that is used by bcuint
  for bicubic interpolation.