This work deals with a geometric inverse source problem. It consists in recovering inclusion in a fixed domain based on boundary measurements. The inverse problem is solved via a shape optimization formulation. Two cost functions are investigated, namely, the least squares fitting, and the Kohn-Vogelius function. In this case, the existence of the shape derivative is given via the first order material derivative of the state problems. Furthermore, using the adjoint approach, the shape gradient of both cost functions is characterized. Then, the stability is investigated by proving the compactness of the Hessian at the critical shape for both considered cases. Finally, based on the gradient method, a steepest descent algorithm is developed, and some numerical experiments for non-parametric shapes are discussed.