A methodology to track bifurcations of periodic orbits in large-scale dissipative systems depending on two parameters is presented. It is based on the application of iterative Newton-Krylov techniques to extended systems. To evaluate the action of the Jacobian it is necessary to integrate variational equations up to second order. It is shown that this is possible by integrating systems of dimension at most four times that of the original equations. In order to check the robustness of the method, the thermal convection of a mixture of two fluids in a rectangular domain has been used as a test problem. Several curves of codimension-one bifurcations, and the boundaries of an Arnold's tongue of rotation number 1/8, have been computed.