I know that the superposition of two Poisson process with rates $\lambda_1$ and $\lambda_2$ is again a Poisson process with rate $\lambda_1+\lambda_2$. Thus this process has interarrival times distributed as an exponential with parameter $\lambda_1+\lambda_2$.
What I am actually doing is numerically sampling the exponential distribution: I collect a bunch of random times sampled from $\lambda e^{-\lambda t}$. The point is that I do it for different values of $\lambda$ and put all together in the same file. Then I plot the distribution. Let's assume the simple case in which I sample an equal number of times for each value of $\lambda$ and that I make $\lambda$ vary say from $10^{-4}$ to $10^3$.
To my intuition the resulting interevent time distribution should be \begin{equation} \int_{\lambda_{min}}^{\lambda_{max}} P(\lambda)\lambda e^{-\lambda t} d\lambda. \end{equation} Given the above hypothesis of uniformly distributed $\lambda$ values I just have $P(\lambda) = const$ and the integral is then \begin{equation} \frac{1}{t}\left(\lambda_{min} e^{-\lambda_{min} t} - \lambda_{max} e^{-\lambda_{max} t}\right) + \frac{1}{t^2}\left(e^{-\lambda_{min} t} - e^{-\lambda_{max} t}\right). \end{equation}
Looking (and plotting) the above function, I expect two different regimes, one in which the interevent time distribution decays as $t^{-1}$ (large times) and one in which decays as $t^{-2}$ (small times). I don't find the second regime in the data.
Question: is it my approach that is wrong or the integral I am doing should actually give me the observed interevent time distribution? If it's correct, why the data don't show the small time regime?
Using a random rate stops it being a Poisson process and, consequently, breaks the independence.
Yes, if you have two (or more) independent Poisson processes, then the superposition of them will be another Poisson process. For homogeneous cases, if they have the intensities $\lambda_1,\dots,\lambda_m$, then the resulting Poisson process is also homogeneous with intensity $\sum_{i=1}^m \lambda_i$.
For the inter-arrival (or inter-occurrence, inter-event etc) times of a single homogeneous Poisson process with intensity $\lambda$, the probability density will be what you have written: $$\lambda e^{-\lambda t}.$$ (Meaning the average inter-arrival time is $1/\lambda$.)
So for the superposition, the probability density of the inter-arrival times is: $$[\sum_{i=1}^m \lambda_i] e^{- [\sum_{i=1}^m \lambda_i]t}.$$
But you are now making $\lambda$ a uniform random variable, say, $U$. This is converting the Poisson process into what is called a doubly stochastic Poisson process or Cox process. Actually, this particular Cox process is called a mixed Poisson process, which the topic of the book by Grandell, because the intensity/mean measure is now a product of a random variable U and the Lebesgue measure (given here by simply $t$ ie length).
The random Poisson measure/parameter breaks the independence between arrivals. It's also discussed in this StatckExchange thread.
But you if consider just the first arrival-time $T_1$, this gives an integral expression: $$ Prob(T_1>t)=E_U[Prob( \text{No arrivals in time }(0,t])]=E_U[e^{-U t} ]= \int_{\lambda_{\min}}^{\lambda_{\max}} p(u) e^{-u t} du , $$ where $p(u)$ is the probability density of $U$. Taking the complement and differentiating with respect to $t$ gives your expression for the probability density of the first arrival-time. Hopefully numerical experiments will support my claims.