Uncertainty propagation for the solution of an integral equation

I have a dataset and I use Maximum Likelihood Estimation to estimate the values of the parameters of a weibull distribution. The MLE theory provides with theoretical Confidence Intervals (asymptotical, or for $n$ samples).

Then, I use the fitted Weibull distribution in an integral equation which is solved numerically :

$Y(t_0) = h(t_0) . int_{0}^{t_0} S(t) dt + S(t_0)$

Where $h$ and $S$ are the hazard function and the survival function of the distribution, and therefore are functions of the parameters.

I would like to propagate uncertainty on the fitted weibull parameters to estimate confidence intervales or quantiles for Y, how could I do that ?
Thanks !