We propose a controllability method for the numerical solution of time-harmonic Maxwells equations in their first-order formulation. By minimizing a quadratic cost functional, which measures the deviation from periodicity, the controllability method determines iteratively a periodic solution in the time domain. At each conjugate gradient iteration, the gradient of the cost functional is simply computed by running any time-dependent simulation code forward and backward for one period, thus leading to a non-intrusive implementation easily integrated into existing software. Moreover, the proposed algorithm automatically inherits the parallelism, scalability, and low memory footprint of the underlying time-domain solver. Since the time-periodic solution obtained by minimization is not necessarily unique, we apply a cheap post-processing filtering procedure which recovers the time-harmonic solution from any minimizer. Finally, we present a series of numerical examples which show that our algorithm greatly speeds up the convergence towards the desired time-harmonic solution when compared to simply running the time-marching code until the time-harmonic regime is eventually reached.