Summary: We investigate numerically models of static spherically symmetric boson-fermion stars in the scalar-tensor theory of gravity with a massive dilaton field. The proper mathematical model of such stars is interpreted as a nonlinear two-parametric eigenvalue problem with an unknown internal boundary. To solve this problem, the Continuous Analogue of Newton Method (CANM) is employed. On each iteration CANM leads to two separate linear boundary value problems with different dimensions inside and outside the star. Along with them a nonlinear algebraic system for the spectral parameters is also solved.
In this way, we determine the behaviour of the basic geometric quantities and functions describing the dilaton field and matter fields building the star. From the physical point of view, the main result is that the structure and properties of such stars depend essentially both on their fermionic and bosonic components.