Solution (source code)

= Solution

Subtract the star's acceleration from that of $M_1$ in an inertial frame. With $\mathbf r_i$ measured from the star,
$$
\ddot{\mathbf r}_1
=-\frac{G(M_*+M_1)}{r_1^3}\mathbf r_1
+GM_2\left[
\frac{\mathbf r_2-\mathbf r_1}{|\mathbf r_2-\mathbf r_1|^3}
-\frac{\mathbf r_2}{r_2^3}
\right].
$$
This has the requested form $\ddot{\mathbf r}_1=\nabla_1(U_1+\mathcal R_1)$ with
$$
\boxed{U_1=\frac{G(M_*+M_1)}{r_1}},
$$
and <disturbing function>
$$
\boxed{\mathcal R_1=GM_2\left[
\frac1{|\mathbf r_2-\mathbf r_1|}
-\frac{\mathbf r_1\mathbin{\cdot}\mathbf r_2}{r_2^3}
\right]}.
$$
The first term is the <direct disturbing function>, the attraction of $M_2$ on $M_1$. The second is the <indirect disturbing function>, which subtracts the acceleration of the star-centered origin by $M_2$.