Presented is a new stochastic algorithm for computing the probability density functions and estimating bifurcations of nonlinear equations of motion whose system parameters are characterized by a Gaussian distribution. Polynomial and Fourier chaos expansions, which are spectral methods, have been used successfully to propagate parametric uncertainties in nonlinear systems. However, it is shown that bifurcations in the time domain are manifested as discontinuities in the stochastic domain, which are problematic for solution with these spectral approaches. Because of this, a new algorithm is introduced based on the stochastic projection method but employing a multivariate B spline. Samples are obtained by choosing nodes on the stochastic axes. These samples are used to build an interpolating function in the stochastic domain. Monte Carlo simulations are then very efficiently performed on this interpolating function to estimate probability density functions of a response. The results from this nonintrusive and non-Galerkin approach are in excellent agreement with Monte Carlo simulations of the governing equations, but at a computational cost 2 orders of magnitude less than a traditional Monte Carlo approach. The probability density functions obtained from the stochastic algorithm provide a rapid estimate of the probability of failure for a nonlinear pitch and plunge airfoil. Nomenclature a h = distance from the airfoil midchord to the elastic axis, positive aft B i;k 1 ;x 1 = the ith B spline of order k 1 for the knot sequence x 1 B i;k 2 ;x 2 = the jth B spline of order k 2 for the knot sequence= the second moment of inertia of the airfoil about the elastic axis, mr 2 b 2 K h = plunge stiffness of the airfoil, m! 2 h K = torsional stiffness of the airfoil about the elastic axis, I ! 2 k 1 , k 2 = order of the B splines on the 1 and 2 axes, respectively m = mass of the airfoil N 1 , N 2 = number of knots on the 1 and 2 axes, respectively S = first moment of inertia of the airfoil about the elastic axis, mx b t = nondimensional time, scaled by V 1 and b V r = reduced velocity, V 1 =! b V 1 = freestream velocity w = Gaussian weight function x = distance from the elastic axis to the center of mass of the airfoil, positive aft x 1 , x 2 = vectors containing the knots on the 1 and 2 axes, respectively = airfoil pitch = airfoil cubic spring constant = airfoil quintic spring constant s = airfoil mass ratio, m= 1 b 2 1 , 2 = zero mean, unit variance Gaussian random variables 1 = freestream density = random basis functions ! r = frequency ratio, ! h =! ! , ! h = uncoupled natural frequencies in pitch and plunge, respectively Superscripts = mean valuẽ = standard deviation = expansion coefficient