If you apply the trapezoidal rule with $N=n, 2n, 4n$ segments, then you get a value
$$
I_N = I_* + c_2 N^{-2} + c_4 N^{-4} + ...
$$
The estimation of the terms of this expansion from two or more numerical evaluations for different $N$ is called Richardson extrapolation or Romberg integration (successively eliminating the terms for $n^{-2k}$).
Now if you compute the error from two integrations for $N=n,2n$ by eliminating the exact value $I_*$, you get
$$
I_n-I_{2n}=c_2n^{-2}(1-2^{-2})+c_4n^{-4}(1-2^{-4})+
$$
so that you get the error estimate
$$
E_n=\frac43(I_n-I_{2n}) = c_2n^{-2} + \frac54 c_4n^{-4}+...
$$
If you now compute
$$
\frac{E_n}{E_{2n}}
=\frac{c_2 + \frac54 c_4n^{-2}+...}{c_22^{-2} + \frac54 c_42^{-4}n^{-2}+...}
=4+\frac{\frac{15}{4}c_4n^{-2}+...}{c_2+\frac5{16}c_4n^{-2}+...}
$$
you can see from the numerical value of this fraction that if it is far away from $4$ that with some high probability there is some error in the implementation.
If you compute a series of these fractions, you should find some middle region of values for $n$ where the quotient is close to $4$. For large values of $n$ one has to add the floating point noise of size $\mu\cdot n$ to the error, it will dominate the error for $n>1/\sqrt[3]\mu\sim 10^5$, where $\mu=2^{-52}\sim 10^{-16}$ is the machine constant of the floating point data type. Thus for large $n$ the error will increase again in a vaguely linear way and the quotient will fluctuate randomly.
If such a middle region does not exist, then the coefficient $c_4$ or a higher one is very large against $c_2$, the function is then either rapidly growing or highly oscillatory and you will need more adapted methods.
See also similar observations in a BDF ODE solver, a fourth order Taylor ODE solver, classical RK4, and one-sided and symmetric and further improved difference quotients.