docs: def jflatten(j):
      """
      Flatten 3_D Jacobian into 2-D.
      """
      nobs, nf, nargs = j.shape
      nrows, ncols = nf * nobs, nargs * nobs
      jflat = np.zeros((nrows, ncols))
      for n in xrange(nobs):
          r, c = n * nf, n * nargs
          jflat[r:(r + nf), c:(c + nargs)] = j[n]
      return jflat