For j=1,…,r, let Pj be a system of generalized prime numbers, NPj the corresponding system of generalized integers, and ζPj(s), s=σ+it, the Beurling zeta-function. In the paper, we consider simultaneous approximation of a collection of analytic functions by shifts (ζP1(s1+iτ),…,ζPr(sr+iτ)). For this, we require that the summatory functions satisfy the bounds MPj(x)−ajx≪Pjxθj with aj>0 and 0⩽θj<1, j=1,…,r. Moreover, we suppose that the mean square of ζPj(s) is bounded for σ∈(σPj,1) with some σPj depending on Pj and θj. Then the main result of the paper asserts that there exists a closed non-empty set FP1,…,Pr of analytic functions such that its functions (f1(s1),…,fr(sr)) are approximated simultaneously by the above shifts of Beurling zeta-functions. It is proved that the set of approximating shifts has a positive lower density (or density with at most countably many exceptions of approximation accuracy). This shows a good joint approximation of Beurling zeta-functions.