In this paper we are concerned with the non-invasive embedding of enriched partition of unity approximations in classical finite element simulations and the efficient solution of the resulting linear systems. The employed embedding is based on the partition of unity approach introduced in Schweitzer and Ziegenhagel (Embedding enriched partition of unity approximations in finite element simulations. In: Griebel M, Schweitzer MA, editors. Meshfree methods for partial differential equations VIII. Lecture notes in science and engineering, Cham, Springer International Publishing; 195-204, 2017) which is applicable to any finite element implementation and thus allows for a stable enrichment of e.g. commercial finite element software to improve the quality of its approximation properties in a non-invasive fashion. The major remaining challenge is the efficient solution of the arising linear systems. To this end, we apply classical subspace correction techniques to design non-invasive efficient multilevel solvers by blending a non-invasive algebraic multigrid method (applied to the finite element components) with a (geometric) multilevel solver (Griebel and Schweitzer in SIAM J Sci Comput 24:377-409, 2002; Schweitzer in Numer Math 118:307-28, 2011) (applied to the enriched embedded components). We present first numerical results in two and three space dimensions which clearly show the (close to) optimal performance of the proposed solver.