We present a new method for calculating the depletion of polymer segments and the depletion induced demixing in solutions of flexible polymers and small globular proteins. The proteins are modelled as hard spheres. The method consists of a mean field calculation of the local polymer volume fraction with the use of a spherical cell-model. Two different types of polymer contribution to the free energy are considered in this calculation. Most of the results presented here are obtained for a low volume fraction expansion of the polymer free energy and an expression similar to the well-known Flory-Huggins free energy. The free energy also contains a squared expression of the gradient in the polymer volume fraction accounting for the steric hindrance suffered by the polymer segments from the spheres and reflecting the entropy loss of the segments in the neighborhood of the surfaces of the spheres. For dilute solutions containing long polymers and either very small or very large spheres, the local polymer volume fractions phi(2)(r) are obtained analytically. Intermediate cases are studied numerically. It is found that phi(2)(r) in the small sphere limit depends on the sphere volume fraction eta, whereas in the large sphere limit phi(2)(r) is eta-independent and tends to the phi(2)(r) near a flat plate. Using the local polymer volume fraction we calculate the macroscopic polymer volume fraction and the free energy. In this calculation we consider besides the cell model a more refined version which also accounts for sphere-sphere correlations. Both cell models yield qualitatively the same results for the spinodals, the quantitative differences between these models being of the order of 5-25 %. The influence of parameters as the polymer chain length, the polymer segment excluded volume, the polymer concentration, the sphere radius on both the depletion profile and the spinodal is discussed. The spinodal, calculated without adaptable parameters, is shown to lie in the neighborhood of the experimental binodal.