Contents 
We consider the continuation of loci of bifurcations of periodic orbits in large scale dissipative systems depending on two parameters obtained after the discretization of PDEs. These kind of computations involve the integration of the variational equations up to order two when NewtonKrylov methods are employed. The actions of the Jacobian of the system, required by the linear solvers, are computed directly, instead of trying to approximate them by finite differences, for accuracy and convergence reasons. As a byproduct the methods allow also the continuation of bifurcations of fixed points. The application to the thermal convection of a binary fluid mixture will be used as an example of application. In our formulation the initial system consists of three partial differential equations for a streamfunction, and the perturbations of the temperature and the concentration of the denser component. This means that a system of twelve PDE's has to be integrated to compute the action of the Jacobian of the system for the bifurcation points. 
