We show that the large N expansion in the multi-trace 1 formal hermitian matrix model is governed by the topological recursion of [24] with extra initial conditions. In terms of a 1d gas of eigenvalues, this model includes - on top of the squared Vandermonde - multilinear interactions of any order between the eigenvalues. In this problem, the initial data (omega(0)(1), omega(0)(2)) of the topological recursion is characterized: for omega(0)(1), by a non-linear, non-local Riemann-Hilbert problem on the discontinuity locus Gamma to determine; for omega(0)(2), by a related but linear, nonlocal Riemann-Hilbert problem on the discontinuity locus Gamma. In combinatorics, this model enumerates discrete surfaces (maps) whose elementary 2-cells can have any topology - omega(0)(1) being the generating series of disks and omega(0)(2) that of cylinders. In particular, by substitution one may consider maps whose elementary 2-cells are themselves maps, for which we propose the name "stuffed maps." In a sense, our results complete the program of the "moment method" initiated in the 90s to compute the formal 1/N expansion in the 1 hermitian matrix model.