Implicit time integration schemes are widely used in computational fluid dynamics to speed-up computations. Indeed, implicit schemes usually allow for less stringent time-step stability constraints than their explicit counterpart. The derivation of an implicit scheme is however a challenging and time-consuming task, increasing substantially with the model equations complexity since this method usually requires a fairly accurate evaluation of the spatial scheme's matrix Jacobian. This work presents a flexible method to overcome the difficulties associated to the computation of the derivatives, based on the forward mode of automatic differentiation using operator overloading. Flexibility and simplicity of the method are illustrated through implicit resolution of various flow models of increasing complexity such as the compressible Euler equations, a two-phase flow model in full equilibrium (Le Martelot et al., 2014) and a symmetric variant (Saurel et al., 2003) of the twophase flow model of (Baer & Nunziato, 1986) dealing with mixtures in total disequilibrium.