I'll give it a start and stick with Fourier methods. The approaches from Carr and Madan (1999) and Fang and Oosterlee (2009) are indeed known to be inaccurate for highly OTM options. I'd suggest to try out one the following three alternatives. The first one seems to be the most relevant one.
Saddlepoint Method
I begin to cite Hirsa (2013):
The saddlepoint method offers considerably better accuracy in pricing out-of-the-money
options than either the FFT or the COS method. However, the algorithm’s accuracy for
at the money and in the money options is somewhat lacking compared to these two other
methods and like the FFT method its solution must be rederived for each different payoff.
Like almost everything in Fourier based option pricing, this idea stems from Carr and Madan (2009). It is based on Madan's et al. (2008) insight that option prices are tail probabilities.
Under the stock measure,
$$\frac{C(S_0,K,T)}{S_0}=\mathbb{E}_0^\mathbb{S}\left[\frac{\max\{S_T-K,0\}}{S_T}\right]=\mathbb{E}_0^\mathbb{S}\left[\left(1-\frac{K}{S_T}\right)^+\right].$$
If $y=\ln\left(\frac{S_T}{K}\right)$ and $\frac{K}{S_T}=e^{-y}$ and $f$ is the density of $\ln\left(\frac{S_T}{K}\right)$ under $\mathbb{S}$, then, using partial integration,
\begin{align*}
\frac{C(S_0,K,T)}{S_0}&=\int_0^\infty (1-e^{-y})f(y)\mathrm{d}y \\
&= 1-\int_0^\infty e^{-y}F(y)\mathrm{d}y \\
&= \int_0^\infty (1-F(y))e^{-y}\mathrm{d}y.
\end{align*}
Given $y$, $1-F(y)$ is just the probability (under $\mathbb{S}$) of the event $\left\{\frac{S_T}{K}>y\right\}$. Furthermore, $e^{-y}$ is just the pdf of an exponential random variable with mean 1. Thus,
\begin{align*}
\frac{C(S_0,K,T)}{S_0}&= P\left[\left\{\ln\left(\frac{S_T}{K}\right)>Y\right\}\right]=P\left[\left\{X-Y>\ln(K)\right\}\right],
\end{align*}
where $X=\ln(S_T)$ under $\mathbb{S}$ and $Y$ an independent exponential random variable.
That's it. The Saddlepoint method is based on the following idea. Let $M$ be the moment generating function of $X-Y$. Then the cumulant generating function is defined via $K(t)=\ln(M(t))$. The Lugannani Rice formula uses CGFs to approximate tail probabilities. That's precisely what we need. Carr and Madan (2009) write down all the equations. Hirsa (2013) also provides a step by step implementation.
I jump straight to the results. This is Table 2.10 from Hirsa.
.
We got Monte Carlo, the COS method, fractional FFT, plain FFT and the saddlepoint method. Parameters are $S_0=100$, $r=3\%$, $\kappa=2$, $\xi=0.5$, $\bar{v}=0.04$, $v_0=0.04$, $\rho=-0.7$ and $T=\frac{1}{2}$. I can state all the equations here if you don't have Hirsa's book. I guess this method is as reasonably accurate as it gets with Fourier methods.
Carr Madan Time Value
Carr and Madan (1999) acknowledge that their FFT method is not flawless and propose an alternative approach that works with time values only. That could improve the pricing of OTM options. The OTM prices are simply
\begin{align}
C^\text{OTM}(S_0,K,T) &= \mathrm{1}_{\{K>S_0\}} \cdot C(S_0,K,T) \\
P^\text{OTM}(S_0,K,T) &= \mathrm{1}_{\{K<S_0\}}\cdot P(S_0,K,T).
\end{align}
Assuming that Fourier transforms exist, you can write
\begin{align*}
C^\text{OTM}(S_0,K,T) &= \frac{1}{\pi}\int_0^\infty \Re\left(e^{-iku} \hat{C}^\text{OTM}(u)\right)\mathrm{d}u, \\
P^\text{OTM}(S_0,K,T) &= \frac{1}{\pi}\int_0^\infty \Re\left( e^{-iku} \hat{P}^\text{OTM}(u)\right)\mathrm{d}u,
\end{align*}
You can combine both equations to $O(S_0,K,T)=C^\text{OTM}(S_0,K,T)+P^\text{OTM}(S_0,K,T)$. For $K>S_0$, $O$ equals the time value of the call option and for $K<S_0$ it is the time value of a put option. The function $O(K)$ is not defined for $K=S_0$. By linearity, $\hat{O}(u)=\hat{C}^\text{OTM}(u) +\hat{P}^\text{OTM}(u)$. Carr and Madan compute that
\begin{align}
\hat{O}(u) &= e^{-rT}\left( \frac{S_0^{i(u-i)}}{i(u-i)}-\frac{S_0^{i(u-i)}e^{(r-q)T}}{iu}-\frac{\varphi_T(u-i)}{u(u-i)} \right), \label{Eqn: Fourier Transform of Time Value}
\end{align}
This is Equation 14 in Carr and Madan (1999), where I add an $S_0\neq1$ and dividends $q\neq0$. The calculation is straight forward, you only need to change the order of integration.
However, the function $O$ is obviously peaks at the money, see below (I use a GBM with $S_0=2$).

This effect gets worse as $T\to0$ because the time value rapidly converges to zero. The function $O$ looks more and more like a Dirac delta function -- and the Fourier transform spreads out further and further (uncertainty principle).
To address this issue Carr and Madan use $\sinh(\alpha(k-\ln(S_0))$ for damping, rather than the standard $e^{\alpha k}$. The Fourier transform of the damped OTM option price is
\begin{align*}
\psi_T(u) &= \int_\mathbb{R} e^{iku} \sinh\left(\alpha \big(k-\ln(S_0)\big)\right)O(k) \mathrm{d}k\\
&= \frac{\hat{O}(u-i\alpha)e^{-\alpha \ln(S_0)}-\hat{O}(u+i\alpha)e^{\alpha \ln(S_0)}}{2},
\end{align*}
which results in
\begin{align*}
O(k) &= \frac{1}{\pi\sinh\left(\alpha \big(k-\ln(S_0)\big)\right)}\int_0^\infty \Re\left(e^{-iku} \psi_T(u)\right) \mathrm{d}u,
\end{align*}
Control Variates
Somewhat related is the idea to subtract the current intrinsic value from the option to obtain the OTM option time value. Then,
\begin{align*}
C(S_0,K,T) = \left(S_0e^{-qT}-Ke^{-rT}\right)^+- \frac{S_0e^{-qT}}{\pi}\int_0^\infty \Re\left( e^{iu\tilde{k}} \frac{\varphi_{X_T}(u-i)-1}{u(u-i)} \right)\mathrm{d}u,
\end{align*}
where $\tilde{k}=\ln\left(\frac{S_0e^{-qT}}{Ke^{-rT}}\right)$.
Cont and Tankov (2004, Chapter 11.1.3) point out that's better to use a smooth function such as the BS call price. Mark Joshi, too, liked using the BS model as a control variate, see also Joshi and Yang (2011). The call option price looks somewhat like this
\begin{align*}
C(S_0,K,T) = C_\text{BS}(S_0,K,T)- \frac{S_0e^{-qT}}{\pi}\int_0^\infty \Re\left( e^{i\tilde{k}u} \frac{\varphi_{X_T}(u-i)-\varphi_{X_T}^\text{BS}(u-i)}{u(u-i)} \right)\mathrm{d}u,
\end{align*}
where $\tilde{k}=\ln\left(\frac{S_0e^{-qT}}{Ke^{-rT}}\right)$. Furthermore, $C_\text{BS}$ denotes the standard Black Scholes (1973) option price and $\varphi_{X_T}^\text{BS}(u)=e^{-\frac{1}{2}\sigma^2 T u(u+i)}$. There's some literature on how to find the optimal volatility for the BS call, see below for the Heston model:
