Let's calculate the derivative straight forward:
\[
\frac{\partial^2 \Phi}{\partial x^2} = \frac{U}{2r_0^2} \frac{\partial^2}{\partial x^2} \left( \alpha x^2 + \beta y^2 + \gamma z^2 \right) = \frac{U}{2r_0^2} \frac{\partial}{\partial x} \left( 2 \alpha x \right) = \frac{U \alpha}{r_0^2}.\]
By analogy we have this relation:
\[ \Delta \Phi = \frac{U}{r_0^2} \left( \alpha + \beta + \gamma \right),\]
so the only way to satisfy $\Delta \Phi = 0$ is to choose coefficients according to
\[ \alpha + \beta + \gamma = 0\]
The relation which connects the wavelength with the frequency is $2\pi /\lambda = \omega / c$. So, the answer is
\[\omega_\text{fr,max}=\frac{2\pi c}{n r_0}\]
The projection of the second Newton's law for the ion onto the $x$-axis is
$$ m \cfrac{d^2x}{dt^2} = ZeE_x = -\cfrac{Zex}{r_0^2}\Big( \alpha U + \tilde{\alpha} V \cos (\omega_{rf} t )\Big)$$
Using the subsistution $\xi = \omega_\text{rf} t /2$ we can obtain
$$ \cfrac{d^2x}{d\xi^2}+ x\Bigg( \frac{4 Ze \alpha U}{m r_0^2 \omega_\text{rf}^2}+\frac{4 Ze \tilde{\alpha} V}{m r_0^2 \omega_\text{rf}^2} \cos 2\xi \Bigg) = 0,$$
so the answer is
\[\xi = \frac{\omega_\text{rf}t}{2}, \quad a_x = \frac{4 Ze \alpha U}{m r_0^2 \omega_\text{rf}^2}, \quad q_x = -\frac{2 Ze \tilde{\alpha} V}{m r_0^2 \omega_\text{rf}^2}\]
By running the Mathieu-stab.py we can obtain a graph that contains some artifacts in the middle of stability region.
Stable motion in the trap means that the motion along each axis is stable. We know that the equations of motion along the $x$- and the $y$-axes are the same but for the $z$-axis we have $a_z = -2a_x$ and $q_z = -2q_z$. This results that motion in the trap is stable when points $(a_x,q_x)$ and $(-2a_x,2q_x)$ are in stability region in the same time.
If the motion along both $x$- and $y$-axes is stable, the ion will pass through the trap, since $E_z=0$. Otherwise the ion will be ejected. So we have to find the intersection of the regions where the motion of the $\rm ^{87}Rm^+$ ion is stable and the motion of the $\rm ^{86}Rb^+$ ion is unstable.
First, let's find the region where solutions of the Mathieu equation for the parameters $(a_x, q_x)$ and $(-a_x, -q_x)$ are stable at the same time. This follows from $a_y = -a_x$ and $q_y = -q_x$.
Now we scale the coordinates $a_x$ and $q_x$ for both isotopes with masses $m_{86} = 86~\text{Da}$ and $m_{87} = 87~\text{Da}$ and draw their stability regions on the same axes.