This work concerns the numerical solution of large-scale systems of nonlinear equations, when derivatives are not available for use, but assuming that all functions defining the problem are continuously differentiable. A hybrid approach is taken, based on a derivative-free iterative method, organized in three phases. The first phase is defined by derivative-free versions of a fixed-point method that employs spectral parameters to define the steplength along the residual direction. The second phase consists on a matrix-free inexact Newton method that employs the GMRES to solve the linear system that computes the search direction. This second phase will only take place if the first one fails to find a better point after a predefined number of reductions in the step size. A third phase, based on directional direct search, should act whenever too many line searches have excessively decreased the steplenght along the inexact-Newton direction. In all stages, the criterion to accept a new point considers a nonmonotone decrease condition upon a merit function. Convergence results are established for two scenarios: for the proposed three phases method and when the last phase is disabled. The numerical performance is assessed through experiments in a set of problems collected from the literature. Both the theoretical and the experimental analysis support the feasibility of the proposed hybrid strategy.