function two_dof_system_implicit_iterative()
    % Define our two-dof system
    k1 = 1;
    k2 = 1;
    k3 = 50;
    F0 = -0.1;
    
    % Set options for the solver
    solverTolerance = 1.0e-10;
    options = optimoptions('fsolve', 'Display', 'iter', 'FunctionTolerance', solverTolerance);

    % Define the linear system of equations
    f=@(u) [F0 - k1*(u(1) - u(2)) - k3*(u(1) - u(2))^3 ;            % Equation 1
               + k1*(u(1) - u(2)) + k3*(u(1) - u(2))^3 - k2*u(2)];  % Equation 2

    % Implicitly SOLVE our system, iteratively.
    u = fsolve(f, [0,0], options);
    
    % Print the output
    disp(' ');
    disp('*******************')
    disp('SOLUTION:');
    disp(['  displacement: u1 = ', num2str(u(1))]);
    disp(['  displacement: u2 = ', num2str(u(2))]);
    disp('*******************')
end
