Summary: Coupled 3D-1D problems arise in many practical applications, in an attempt to reduce the computational burden in simulations where cylindrical inclusions with a small section are embedded in a much larger domain. Nonetheless the resolution of such problems can be non trivial, both from a mathematical and a geometrical standpoint. Indeed 3D-1D coupling requires to operate in non standard function spaces, and, also, simulation geometries can be complex for the presence of multiple intersecting domains. Recently, a PDE-constrained optimization based formulation has been proposed for such problems, proving a well posed mathematical formulation and allowing for the use of non conforming meshes for the discrete problem. Here an unconstrained optimization formulation of the problem is derived and an efficient gradient based solver is proposed for such formulation, along with a suitable preconditioner to speed up the iterative solver. Some numerical tests on quite complex configurations are discussed to show the viability of the method .