We present an analytical solution based on the series expansion of Papkovich-Boussinesq's displacement potentials to derive the elastic solution of a spherical inclusion containing an eigenstrain or pressurized cavity in an elastic half-space. The inclusion–host interface can be treated as perfectly bonded or allowing different amount of interface sliding as a function of shear stress. The analytical solution allows systematic investigation on some key parameters that control the elastic stress field, such as inclusion depth, elastic moduli and the amount of interface sliding. The model is applied to study the distributions of stresses and displacements within and around the inclusion. Stress trajectories and slip lines are computed around a pressurized cavity based on the analytical solution to study potential fracture modes and patterns. The amount of inclusion pressure relaxation due to the free surface is also systematically investigated as a function of inclusion depth and shear modulus ratio between the host and inclusion. A MATLAB code is provided that allows one to directly apply the analytical solution to natural systems given any elastic parameters. The code is benchmarked with Mindlin's solution for loaded homogeneous inclusion in proximity to a free surface and 3-D finite element simulations for heterogeneous inclusion. This model may contribute to the study of the mechanical properties of natural systems at various scales, from km size magma chamber to mm-μm size mineral inclusions sealed in thin section.