We develop a new, efficient, and accurate method of obtaining numerical solutions for soliton bag models. By introducing additional differential equations to describe the integral constraints and eigenvalue conditions, and after a careful analysis of the asymptotic form of the solutions, we convert the equations of the soliton bag model in the mean-field approximation into a two-point boundary-value system, which is readily solved by standard methods. To illustrate the techniques we find solutions for a chiral model, for cases of both massless and massive pions, for a wide range of parameters, and compare our results with previous calculations.