This is not a question about MATLAB, but a question about understanding numerical integration, and perhaps a bit of statistics. Not only that, you need to appreciate that the result that trapz produces is only an approximation, because if you knew the EXACT values at each point, trapz will return only an approximation to the "true" integral of that function. In fact, the true integral could be almost anything, since we don't know what the function does between those points.
So at best, all you can ask is the uncertainty on the result that trapz would have generated, If the data had error bars of zero width at each point. Again, that still has nothing to do with MATLAB, because we still have not yet defined the meaning of your question.
So lets assume that we have a set of points, as (x,y) pairs. Assume that x is a monotonic set of values, so x(i+1) > x(i). At each point, y(i) is given, but assume a set of error bars around y(i) exist and are known.
Now, let me consider what an error bar means. Does that say that the true function is uniformly ANY value within the limits of those error bars? I'll assume that is true, but you may prefer to assume the function is normally distributed, with say 99% of the values lay inside the error bars. That would change the analysis in a somewhat minor way.
Regardless, the important thing to remember is that integration is a LINEAR operator! What does that mean? That tells us
int( F(x) + G(x) ) = int(F(x)) + int(G(x))
Here int can be used as a substitute for trapz, or any such tool.
So assume the function is given as a set of limits, lower and upper bounds on the function at each point. So now we have a vector x, and vectors ylo and yhi. The analysis is simpler if the vector x is uniformly spaced, since then trapezoidal integration is way simpler to write.
So next, we can write the problem as a SUM of two integrals. That is, break down the problem as the integral:
int(ylo(x)) + int(e(x))
where e(x) is a strictly positve, but random function that varies uniformly beween the limits 0 and yhi(x)-ylo(x).
Now we can try to discuss the distribution of the result. How do those limits vary? Are the width of those error bars uniform, so all the same size? If so, then yhi(x)-ylo(x) is a constant. In fact, now the central limit theorem starts to come into the problem. We can now compute the mean and variance of the result.
All of this is pure speculation, since I've now made a LOT of assumptinos, only some of which are probably valid.
However, I would suggest that you talk to someone with an understanding of both numerical integration and at least some knowledge of statistics. While I am an excellent person in that respect, I don't do consulting. (So please don't send me e-mail asking for me to help you directly. I won't.)