Here's the straight answer you're looking for:
You can't calculate \$P_{-1dB}\$ from \$IP_2\$.
Why? Because \$IP_2\$ doesn't contain any information about the cubic term of the non-linearity (\$a_3\$), as it only depends on the linear (\$a_1\$) and quadratic terms (\$a_2\$). Recall that:
$$
s_o(t)=a_1s_i(t)+a_2s_i^2(t)+a_3s_i^3(t)
$$
For a single-tone input \$s_i(t)=A_0cos(\omega_0t)\$ where \$\omega_0=2{\pi}f_0\$, this is what you get at the output:
$$
\begin{align}
s_o(t)&=s_{o,\ DC}(t)+s_{o,\ f_0}(t)+s_{o,\ 2f_0}(t)+s_{o,\ 3f_0}(t)\\
\text{where:}\\
s_{o,\ DC}(t)&=\frac{1}{2}a_2A_0^2\\
s_{o,\ f_0}(t)&=\left(1+\frac{3}{4}\frac{a_3}{a_1}A_0\right)a_1A_0cos(\omega_0t)\\
s_{o,\ 2f_0}(t)&= \frac{1}{2}a_2A_0^2cos(2\omega_0t)\\
s_{o,\ 3f_0}(t)&=\frac{1}{4}a_3A_0^3cos(3\omega_0t)
\end{align}
$$
The gain compression is a consequence of the cubic term (more generally speaking: of the odd-order terms). That's why you can calculate \$P_{-1dB}\$ from \$IP_3\$: because it contains information about \$a_3\$ (and about \$a_1\$, as well), but not from \$IP_2\$:
$$
\begin{align}
IP_{i,2}&=\frac{a_1}{a_2}\\
IP_{i,3}&=\sqrt{\frac{4}{3}\left\lvert\frac{a_1}{a_3}\right\rvert}\\
\end{align}
$$
Also note that your Matlab code doesn't even consider the cubic term:
% Parameters
a1=5; a2=4;