Abstract.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. 1. Introduction. The study of dynamical systems involves, in addition to pure numerical simulations, the computation of invariant manifolds (fixed points, periodic orbits, invariant tori, etc.) and the connections among them (homo-and heteroclinic orbits and heteroclinic chains), the investigation of the stability of these objects, and the examination of their bifurcations when the parameters present in the system are varied. These essential tools help to understand the full dynamics of the system and its dependence on the parameters. The existence of robust continuation and bifurcation packages such as AUTO [10], CONTENT [22], MATCONT [7], etc., allows many of these calculations to be almost routinely performed for moderate-dimensional systems of ordinary differential equations (ODEs). These packages include the continuation of codimension-one bifurcations of fixed points, and even of periodic orbits, and some also include the detection and analysis of codimension-two points. They implement direct solvers for the linear systems involved in the computations, and find the full spectrum when solving eigenvalue problems to study of the stability of the invariant objects.The theory of the extended systems used in these packages to follow bifurcations of fixed points is well developed and can be found in, among others, [42,26,18,53,29,5,16]. The bordered systems for periodic orbits, based on boundary value problems, are analyzed in [11]. In this latter case piecewise collocation in time is used instead of shooting methods.The difficulties in the application of these methodologies to high-dimensional systems, obtained in most cases by discretizing systems of partial differential equations (PDEs), come