Release kinetics for heterogeneous sphere ensembles with a dissolved drug, i.e., initial drug loading below or equal to the drug solubility in the matrix, in a finite external medium was modeled with consideration of heterogeneity among and within spheres. Numerical solutions were obtained using the finite element method for sphere ensemble with normal or log-normal distribution of particle size or initial drug loading among spheres. Exact series solutions were derived for ensembles with various initial loading distributions within spheres, namely linear, quadratic, sigmoidal and uniform distribution, using their mean or average radii. Simplified solutions retaining only one term of the series for non-uniform distributions and three terms for uniform distribution were suggested because of their good approximation to the exact solution. The results of finite element analysis showed that the release rate of an ensemble decreased with increasing standard deviation of particle size. Using weight-average radii in the exact solution gave a prediction of release profile closer to that from the actual size distribution than using mean radii. The three non-uniform loading patterns within spheres all showed reduced initial burst and release rate, leading to more steady release rates than uniform loading, among which the sigmoidal distribution offered the best near-zero order release. Non-uniform initial loading among spheres seemed to have insignificant influence on the release profiles. The volume ratio of liquid to a sphere ensemble played an important role in release kinetics. The derived analytical solutions are applicable to multiple spheres or a single sphere in a finite medium or in a perfect sink.