Variability in multiple independent input parameters makes it difficult to estimate the resultant variability in a system's overall response. The propagation of errors (PE) and Monte Carlo (MC) techniques are two major methods to predict the variability of a system. However, the formalism of PE can lead to an inaccurate estimate for systems that have parameters varying over a wide range. For the latter, the results give a direct estimate of the variance of the response, but for complex systems with many parameters, the number of trials necessary to yield an accurate estimate can be sizeable to the point the technique becomes impractical. The effectiveness of a designed experiment (orthogonal array) methodology, as employed in Taguchi tolerance design (TD) method to estimate variability in complex systems is studied. We use a linear elastic three-point bending beam model and a nonlinear extended finite elements crack growth model to test and compare the PE and MC methods with the TD method. Our focus is on assessing their effectiveness under those specific circumstances. Results from an MC estimate, using 10,000 trials, serve as a referent to evaluate the result in both cases. We find that the PE method works suboptimal for a coefficient of variation above 5% in the input variables. In addition, we find that the TD method works very well with moderately sized trials of designed experiments for both models. Our results demonstrate how the variability estimation methods perform in the deterministic domain of numerical simulations and can assist in designing physical tests by providing a guideline performance measure.