“> Is it really a big deal to solve ∇×u=f in 2D using finite elements ?”
You really want to base a formulation on some type of variational form, which exhibits useful stationarity properties. I suggest you have a look at a finite element reference to see why this is. Silvester and Ferrari’s “Finite Elements for Electrical Engineers” is a good place to start.
In brief, you generally want to start with the functional
F(\mathbf{A}) = \frac{1}{2}\iint\limits_{\Omega}\left\lbrace\nabla\times\mathbf{A}\cdot\nabla\times\mathbf{A}+\sigma\mathbf{A}\cdot\frac{\partial\mathbf{A}}{\partial t}\right\rbrace\, dS + B(\mathbf{A},\mathbf{B}^\mathrm{inc})_{\partial\Omega},
where B() is a boundary operator over the perimeter \partial\Omega where you treat the applied field similar to how you would treat a scattering problem.
If you really want to treat \nabla\times\mathbf{u}=\mathbf{f} directly, you could always form the least-squares residual,
R=\vert\!\vert\nabla\times\mathbf{u-f}\vert\!\vert^2
and find the solutions that minimize R. Note that this changes how you treat the boundary conditions, though.