If your random variables X and Y are independent, Cupid might be helpful. It defines a general convolution RV for any continuous independent RVs. If your distributions are already defined in cupid, the code would look something like:
X = Exponential(1/10); % Illustration with example distributions
Y = Beta(5,2.3);
Z = Convolution(X,Y);
a = 2;
b = 4;
probInterval = Z.CDF(b) - Z.CDF(a);
If either/both of your distributions are not already in Cupid, it is pretty easy to add new distributions if you have the formula for either the PDF or the CDF.
Whether this will be fast and accurate enough for you, I cannot say.
