We study semi-linear elliptic PDEs with polynomial non-linearity and provide a probabilistic representation of their solution using branching diffusion processes. When the non-linearity involves the unknown function but not its derivatives, we extend previous results in the literature by showing that our probabilistic representation provides a solution to the PDE without assuming its existence. In the general case, we derive a new representation of the solution by using marked branching diffusion processes and automatic differentiation formulas to account for the non-linear gradient term. In both cases, we develop new theoretical tools to provide explicit sufficient conditions under which our probabilistic representations hold. As an application, we consider several examples including multi-dimensional semi-linear elliptic PDEs and estimate their solution by using the Monte Carlo method.