We compute families of spherically symmetric neutron-star models in two-derivative scalar-tensor theories of gravity with a massive scalar field. The numerical approach we present allows us to compute the resulting spacetimes out to infinite radius using a relaxation algorithm on a compactified grid. We discuss the structure of the weakly and strongly scalarized branches of neutron-star models thus obtained and their dependence on the linear and quadratic coupling parameters α0, β0 between the scalar and tensor sectors of the theory, as well as the scalar mass μ. For highly negative values of β0, we encounter configurations resembling a “gravitational atom”, consisting of a highly compact baryon star surrounded by a scalar cloud. A stability analysis based on binding-energy calculations suggests that these configurations are unstable and we expect them to migrate to models with radially decreasing baryon density and scalar field strength.