Suppose that your predictive density for $X$ is uniform on $[0,1]$, while your predictive density for $Y$ is standard lognormal, and $f(x,y)=x^y$.
Then you can draw (a lot of) random numbers from each predictive density, plug them into $f$ and get a predictive distribution for $Z$, from which you can deduce expected errors (any error measure you want) and prediction intervals. In R:
nn <- 10000
xx <- runif(nn)
yy <- rlnorm(nn)
Note that for this to work, you will need to specify the full predictive densities - prediction intervals alone won't help you. That is inherent in the problem. Any prediction interval for (say) $X$ can be consistent with many predictive densities for $X$, and with many different resulting distributions for $Z$.