code: private static double[] newtonSolve(Function<double[], double[]> function, double[] initial, double tolerance) {
          RealVector dx = new ArrayRealVector(initial.length);
          dx.set(tolerance + 1);
          int iterations = 0;
          int d = initial.length;
          double[] values = Arrays.copyOf(initial, initial.length);
  
          while (dx.getNorm() > tolerance) {
              double[] fx = function.apply(values);
              Array2DRowRealMatrix df = new Array2DRowRealMatrix(fx.length, d);
              ArrayRealVector fxVector = new ArrayRealVector(fx);
              for (int i = 0; i < d; i++) {
                  double originalValue = values[i];
                  values[i] += H;
                  double[] fxi = function.apply(values);
                  values[i] = originalValue;
                  ArrayRealVector fxiVector = new ArrayRealVector(fxi);
                  RealVector result = fxiVector.subtract(fxVector);
                  result = result.mapDivide(H);
                  df.setColumn(i, result.toArray());
              }
              dx = new RRQRDecomposition(df).getSolver().solve(fxVector.mapMultiply(-1));
              // df has size = initial, and fx has size equal to whatever that function produces.
              for (int i = 0; i < values.length; i++) {
                  values[i] += dx.getEntry(i);
              }
              iterations++;
              if (iterations % 100 == 0) {
                  tolerance *= 10;
              }
          }
          return values;
      }
docs: Finding the best fit using least-squares method for an equation system
  
  @param function Equation system to find fit for. Input: Parameters, Output: Residuals.
  @param initial Initial 'guess' for function parameters
  @param tolerance How much the function parameters may change before a solution is accepted
  @return The parameters to the function that causes the least residuals