In this paper we propose an extension of the generalized Lagrangian multiplier method (GLM) of Munz et al. [52,30], which was originally conceived for the numerical solution of the Maxwell and MHD equations with divergence-type involutions, to the case of hyperbolic PDE systems with curl-type involutions.The key idea here is to solve an augmented PDE system, in which curl errors propagate away via a Maxwelltype evolution system. The new approach is first presented on a simple model problem, in order to explain the basic ideas. Subsequently, we apply it to a strongly hyperbolic first order reduction of the CCZ4 formulation (FO-CCZ4) of the Einstein field equations of general relativity, which is endowed with 11 curl constraints. Several numerical examples, including the long-time evolution of a stable neutron star in anti-Cowling approximation, are presented in order to show the obtained improvements with respect to the standard formulation without special treatment of the curl involution constraints.The main advantages of the proposed GLM approach are its complete independence of the underlying numerical scheme and grid topology and its easy implementation into existing computer codes. However, this flexibility comes at the price of needing to add for each curl involution one additional 3 vector plus another scalar in the augmented system for homogeneous curl constraints, and even two additional scalars for non-homogeneous curl involutions. For the FO-CCZ4 system with 11 homogeneous curl involutions, this means that additional 44 evolution quantities need to be added.Keywords: generalized Lagrangian multiplier approach (GLM), hyperbolic PDE systems with curl involutions, Einstein field equations with matter source terms, first order reduction of the CCZ4 system (FO-CCZ4), stable neutron star in anti-Cowling approximation