In the paper a probabilistic methodology based upon an approximation of the stability region boundary is presented. The proposed approach, taking into account the randomness of the main variables and parameters, such as fault clearing time, the reclosing time and fault location, is able to determine the actual transient stability margins. The analyzed problem has been formalized as a particular free boundary-optimization problem whose resolution is addressed by the adoption of an implicit numerical method, specifically Heun's method. Furthermore, a novel probability stability index is defined. A numerical application validates the aforementioned methodology in the last section of the paper.