pdfs are continuous functions, so the closer spaced your x points are, the closer you get to the expected answer. In the following code the smaller you make the array spacing delx, the closer you get to the expected answer (although for delx = .0001 it takes awhile).
n1 = 3; n2 = 5;
delx = .001;
xmax = -ceil(log(1e-10)*n1)
x = 0:delx:xmax;
a1 = exp(-x/n1)/n1;
a2 = exp(-x/n2)/n2;
uh = conv(a1,a2)*delx;
x3 = 0:delx:2*xmax;
I don't have the exppdf function so I used the equivalent.
Trapz, primitive as it is, shows the trend, but it would be much preferable if Mathworks supported something better in basic Matlab like integration by cubic spline.
 Convolution of exponential pdfs
In the case of the convolution of n exponential pdfs, all of whose decay constants are unique (no repeats): let 'a' be the vector of decay constants. The kth pdf is
The convolution of all n pdfs is the sum over k of
c(k) = a(k)^(n-1) / [ (a(k)-a(1))*(a(k)-a(2)) ...
In the denominator the factor (a(k) - a(k)) = 0 is excluded, so there are n-1 factors in all.
 Convolution by fft
For a convolution of n general pdfs, let x be a row array of N equally spaced points with spacing delx, where the range of x is wide enough that all pdfs die down to very small values at the upper and lower ends of the x array. Let M be an (n x N) matrix of the pdfs stacked on top of each other. Then the convolution is
y = real(ifft(prod(fft(M,,2))))*delx^(n-1);
In other words take the fft of each pdf down the rows, multiply them all together down the columns and transform back. For n = 10 and N = 2e6 points this takes less than 2 seconds on my PC and it is basically linear in N (as long as N has lots of divisors).