% The system of non-linear equations
f1 = @(x1, x2) x1.^3 + x2.^3 - 8 * x1 * x2;
f2 = @(x1, x2) x1 * log(x2) - x2 * log(x1);

% The target function
Fi = @(x1, x2) f1(x1, x2).^2 + f2(x1, x2).^2;

% Accuracy
epsilon = 0.001;
% Maximum iterations
maxits = 100;
% Number of unknowns
unknownsNumber = 2;

% Step for the grid
shiftingStep = 1;
% Grid  X Y
grid = [1 1   % start
        3 3]; % end

% What start point must be shifting
mustShiftGridX = false;

% Stabilization for first shift condition
grid(1, 2) = grid(1, 2) - shiftingStep;

fprintf(' x1     x2     x1:beg x2:beg\n');

% Break if grid so small
 while grid(1, 1) < grid(2, 1) || grid(1, 2) < grid(2, 2)

      % Shifting for each dimension
      if mustShiftGridX == true
           grid(1, 1) = grid(1, 1) + shiftingStep;
      else
           grid(1, 2) = grid(1, 2) + shiftingStep;
      end
      mustShiftGridX = ~mustShiftGridX;

      % The start approximation
      XYaprx = [grid(1, 1) grid(1, 2)];
      detL = 0.5;

      X = zeros(3, 2); % The temp value like a lyambda
      L = zeros(1, 3); % Lyambda
      a = zeros(1, 3); % Result of system equations
      F = zeros(1, 3); % Result of Fi functions
      e = [1 0
           0 1];

      % Break if count of iterations more than maximum
      for i = 1:maxits
        A = Fi(XYaprx(1), XYaprx(2));
        for j = 1:unknownsNumber

            % Lyambdas and lyambda results calculations
            L(2) = detL;
            X(1, :) = XYaprx;
            X(2, :) = XYaprx + L(2) * e(j, :);
            F(1) = Fi(X(1, 1), X(1, 2));
            F(2) = Fi(X(2, 1), X(2, 2));
            if F(2) < F(1)
                L(3) = L(1) + 2 * detL;
            else
                L(3) = -detL;
            end
            X(3, :) = XYaprx + L(3) * e(j, :);
            F(3) = Fi(X(3, 1), X(3, 2));

            % Finding a's
            a(1) = ( (F(2) - F(1)) * L(3) - (F(3) - F(1)) * L(2) ) / ( L(2)^2 * L(3) - L(2) * L(3) );

            a(2) = ( (F(3) - F(1)) * L(2)^2 - (F(2) - F(1)) * L(3)^2 ) / ( L(2)^2 * L(3) - L(2) * L(3) );

            a(3) = F(1);

            Lmin = -a(2) / (2 * a(1));
            % Make a step
            XYaprx = XYaprx + Lmin * e(j, :);
        end
            B = Fi(XYaprx(1), XYaprx(2));
            if abs(A - B) <= epsilon
                break;
            else
                % Decrease the local step
                detL = detL / 2;
                j = 1;
            end
      end
    fprintf(' %.4f %.4f %.4f %.4f \n', XYaprx(1), XYaprx(2), grid(1, 1), grid(1, 2));
 end