Incorporating probabilistic terms in mathematical models is crucial for capturing and quantifying uncertainties in real-world systems, especially when the solution is not unique or exhibits sudden qualitative changes as parameters vary. However, stochastic models typically require large computational resources to produce meaningful statistics. In this work, we leverage the Polynomial Chaos (PC) expansion to propose a systematic approach for bifurcation detection in parametric systems of equations. We show that the method, exploiting a perturbed version of the deterministic model, avoids repeated costly simulations across multiple parameter values and requires no prior information for initializing numerical solvers, while still providing accurate characterization of the bifurcation branches. We argue that the PC solutions of the perturbed model not only provide access to statistical information about the deterministic branches, but also approximate simultaneously their behavior in a single solver call. Finally, we validate our claims by means of numerical tests on the pitchfork bifurcation, examining both its normal form and a classical realization in fluid-dynamics PDEs, namely the Coandà† effect.
