There are no nontrivial examples with $n\ge3$ beyond what I mentioned in my comment above, namely, either a conformal connected sum of a compact space form $(M_1,g_1)$ with the standard round $n$-sphere with the connected sum $M_1\# M_2$ being homothetic to $(M_1,g_1)$ or the case where both $M_1$ and $M_2$ are diffeomorphic to a (possibly exotic) $n$-sphere with homothetic metrics $g_1$ and $g_2$ and the resulting connected sum being again homothetic to $(M_1,g_1)$. (Note that there are many inequivalent Einstein metrics on manifolds homeomorphic to odd-dimensional spheres.)
To begin, consider the following question: Given an Einstein $n$-manifold $(M^n,g)$ (that is connected but not necessarily complete), with Einstein constant $\lambda$ (i.e., $\mathrm{Ric}(g)= (n{-}1)\lambda g$), what are the possibilities for Einstein metrics on $M$ that are conformally equivalent to $g$?
By the well-known formula for the Ricci tensor of $\tilde g = u^{-2}\,g$, where $u>0$ is a function on $M$, $$ \mathrm{Ric}(\tilde g) = \mathrm{Ric}(g) + (n{-}2)\,u^{-1}\,\nabla(\mathrm{d}u) - \bigl(u^{-1}\,\Delta u + (n{-}1)\, u^{-2}\,|\mathrm{d}u|^2\bigr)\, g\,, $$ it follows that $\tilde g$ is also Einstein if and only if $u$ is a non-vanishing function on $M$ that satisfies the equation $\nabla(\mathrm{d}u) = v\,g$ for some function $v$. Differentiating this equation, one finds that $v+\lambda u$ must be locally constant. Thus, by the connectedness of $M$, there exists a constant $c$ such that $\nabla(\mathrm{d}u) = (c{-}\lambda u)\,g$. Consequently, $\Delta u = n(\lambda u{-}c)$, which will be useful below. Moreover, one sees that the Einstein constant of $\tilde g$ is $\tilde\lambda = 2c\,u - \lambda\,u^2 - |\mathrm{d}u|^2$.
(Also, one sees that the vector field $\nabla u$ must annihilate the Weyl tensor of $g$, i.e., $W(\nabla u, X,Y,Z) = 0$ for all vector fields $X,Y,Z$ on $M$. However, that will play no role in this argument.)
For example, if $(M^n,g)=\bigl(\mathbb{R}^n,|\mathrm{d}x|^2\bigr)$ is the standard flat metric, then $\tilde g = u^{-2}\,g$ is Einstein (on the open set where $u$ is nonzero) if and only if $u = p\,|x|^2 + 2\,q{\cdot}x + r$ for some constants $p,r\in\mathbb{R}$ and vector $q\in\mathbb{R}^n$, not all zero.
This shows how to essentially 'linearize' the problem of describing the conformal multiples of an Einstein metric that are, themselves, Einstein: On the bundle $E = \mathbb{R}\oplus\mathbb{R}\oplus T^*\!M$ over $M$, consider the connection $D$ that, for any vector field $X$ and any section $(c,u,\alpha)$ of $E$, satisfies $$ D_X\begin{pmatrix}c, & u, & \alpha\end{pmatrix} = \begin{pmatrix}\mathrm{d}c(X),& \mathrm{d}u(X)-\alpha(X),& \nabla_X\alpha - (c-\lambda\,u) X^\flat\end{pmatrix} $$ where $X^\flat$ is the $1$-form that satisfies $X^\flat(Y) = g(X,Y)$ for all vector fields $Y$ on $M$. Then we can rephrase the above equation on $u$ (and $c$) by saying that $u^{-2}\,g$ is Einstein if and only if there is a $c$ such that $(c,u,\mathrm{d}u)$ is a $D$-parallel section of $E$.
Now, it turns out that, except for the metrics of constant sectional curvature on the $n$-sphere, for a compact, connected Einstein manifold $(M,g)$, the space of global $D$-parallel sections of $E$ consists only of the sections of the form $(\lambda u,\,u,\,0)$ where $u$ is constant. (Such sections are obviously $D$-parallel.)
To see this, first note that if $(M^n,g)$ is compact, connected and Einstein with a non-positive Einstein constant, then the only global $D$-parallel sections of $E$ are those with $u$ constant. This is because the defining equation implies $\Delta u = n(\lambda u {-} c)$. When $\lambda=0$, this clearly implies $c=0$ and that $u$ be constant. When $\lambda<0$, this implies that $\lambda u {-} c$ is an eigenvalue of $\Delta$ with negative eigenvalue $n\lambda$ and hence must be $0$.
Meanwhile, if one supposes that $\lambda>0$, the equation $\nabla(\mathrm{d}u) = (c{-}\lambda u)\,g$ implies that, if $p$ is a critical point of $u$, then the Hessian of $u$ at $p$ is $\bigl(c{-}\lambda u(p)\bigr)$ times the metric at $p$. If $c{-}\lambda u(p)=0$, then the $D$-parallel section $(c,u,\mathrm{d}u)$ agrees with the $D$-parallel section $(\lambda\,u(p), u(p),0)$ at $p$ and hence must equal it, i.e., $u$ must be constant. On the other hand, if $c{-}\lambda u(p)\not=0$, then $p$ is a non-degenerate critical point of $u$ that is either a maximum or a minimum. Thus, on a connected, compact Einstein manifold $(M^n,g)$ the only critical points of a nonconstant $u$ such that $(c,u,\mathrm{d}u)$ is $D$-parallel for some $c$ are either strict maxima or minima. By the Mountain Pass Lemma, it follows that $u$ can only have one maximum and one minimum, which implies that $M$ is homeomorphic to a $n$-sphere.
Suppose that such a non-constant $u$ exists on $(M^n,g)$. By adding a constant to $u$, we can suppose that $u$ is strictly positive on $M$ and by scaling, we can assume that $\lambda=1$. Now by using the structure equations, one can prove that, on a punctured neighborhood of a maximum $p\in M$, $g$ can be written in the 'sine-cone' form $g = \mathrm{d}r^2 + (\sin r)^2\, h$, where $r$ is the distance from $p$, $h$ is an Einstein metric on $S^{n-1}$ with Einstein constant $1$, and $u = a + b \cos r$ for some constants $a$ and $b$. However, unless $h$ is conformally flat, the above sine-cone metric cannot be smooth at $r=0$. Thus, $h$ is conformally flat, which implies that $g$ is also conformally flat. Hence, the only possibility for $(M,g)$ is the round $n$-sphere, as was to be shown.
Now, for any connected Einstein $(M^n,g)$, the sheaf of $D$-parallel sections of $E$ is well-behaved: There is an integer $k$ satisfying $1\le k\le n{+}2$ such that, for any 1-connected open set $U\subset M$, the dimension of the vector space $\mathcal{P}(U)\subset\Gamma(U,E)$ consisting of $D$-parallel sections of $E$ over $U$ is $k$. This follows from the well-known result of DeTurck and Kazdan that $g$, being Einstein, is real-analytic in harmonic coordinates: Since the linear connection $D$ is perforce real-analytic, any germ of a $D$-parallel section defined on an open neighborhood of $p\in M$ can be real-analytically continued as a $D$-parallel section along any real-analytic curve that starts at $p$.
One more useful fact that is easy to prove: If $(c,u,\mathrm{d}u)$ is a $D$-parallel section of $E$ on a connected domain $U\subset M$ that is not identically zero, then the zero locus of $u$ is either empty; a nonempty smooth hypersurface, in which case $\tilde\lambda < 0$; or consists of isolated points that are non-degenerate maxima or minima, in which case, $\tilde\lambda = 0$.
With all this understood, let's look at what might reasonably described as a 'conformal connected sum' of compact Einstein manifolds in dimension $n\ge 3$. I'll take it to be this: Three compact Einstein $n$-manifolds $(M_i,g_i)$ for $i=0, 1, 2$, with respective Einstein constants $\lambda_i$, together with smoothly embbeded compact submanifolds $N_i\subset M_i$ where $N_0$ is diffeomorphic to $[0,1]\times S^{n-1}$, and $N_i$ is diffeomorphic to the $n$-ball $B^n$ for $i=1,2$, and, finally, where $M_0{\setminus}N_0$ is conformally diffeomorphic (as Riemannian manfolds) to the disjoint union of $M_1{\setminus}N_1$ and $M_2{\setminus}N_2$. (Without lost of generality, one can take this diffeomorphism to be the identity map, so assume this.)
Given these hypotheses, there exist smooth functions $u_i: M_i{\setminus}N_i\to (0,\infty)$ for $i=1,2$, such that $g_0 = {u_i}^{-2} g_i$ on $M_i{\setminus}N_i\subset M_0{\setminus}N_0$. Thus, there exist constants $c_i$ such that $(c_i,u_i,\mathrm{d}u_i)$ is a $D_i$-parallel section of $E_i$ on $M_i{\setminus}N_i$. Since $u_i$ is defined on the complement of the smoothly embedded ball $N_i\subset M_i$ in $M_i$, the 1-connectedness of $N_i$ implies that this $D_i$-parallel section extends uniquely as a $D_i$-parallel section of $E_i$ on all of $M_i$. In particular, $u_i$ extends to all of $M_i$ satisfying $\nabla(\mathrm{d}u_i)=(c_i{-}\lambda_i\,u_i)g_i$ for $i=1,2$.
Consequently, for $i=1$ or $2$, $u_i$ must be constant unless $(M_i,g_i)$ is the round $n$-sphere.
If both $u_1$ and $u_2$ are constant, then by scaling $g_1$ and $g_2$ by constants, we can assume that $g_0 = g_i$ on $M_i{\setminus}N_i\subset M_0{\setminus}N_0$ for $i=1$ and $2$, and hence $\lambda_0=\lambda_1=\lambda_2$. Moreover, because $N_i$ is simply-connected and $g_0$ and $g_i$ are real-analytic, it is not hard to show that the isometric inclusion $M_i{\setminus}N_i\subset M_0{\setminus}N_0$ extends to an isometry $\iota_i:(M_i,g_i)\to (M_0,g_0)$, and the compactness of $M_i$ and $M_0$ then implies that $\iota_i$ is a diffeomorphism, mapping $N_i$ diffeomorphically onto $M_0{\setminus}M_i = N_0\cup (M_{3-i}{\setminus}N_{3-i})$. Thus, $M_0$, $M_1$, and $M_2$ are all isometric and homeomorphic to an $n$-sphere. While $(M_0,g_0)$ is indeed a 'conformal connected sum' of $(M_1,g_1)$ and $(M_2,g_2)$ in the above sense, these ought to be regarded as trivial cases. As remarked at the beginning, though, there are many, many such examples that are not conformally flat.
If, say, $u_1$ is constant and $u_2$ is not, we can still reduce to the case that $g_0=g_1$ on $M_1{\setminus}N_1\subset M_0{\setminus}N_0$, and hence $\lambda_0=\lambda_1$. Again, because $N_1$ is simply-connected and $g_0$ and $g_1$ are real-analytic, it follows that the isometric inclusion $M_1{\setminus}N_1\subset M_0{\setminus}N_0$ extends to an isometry $\iota_1:(M_1,g_1)\to (M_0,g_0)$, and the compactness of $M_1$ and $M_0$ then implies that $\iota_1$ is a diffeomorphism, mapping $N_1$ diffeomorphically onto $M_0{\setminus}M_1 = N_0\cup (M_2{\setminus}N_2)$. Meanwhile, $(M_2,g_2)$ must be isometric to the round $n$-sphere, which is conformally flat, implying that $(M_0,g_0)$ is conformally flat. Of course, examples such as this do exist with all of the $(M_i,g_i)$ being conformally flat, and $u_2$ non-constant, in which one attaches a conformally flat bubble to a compact space form, but $(M_0,g_0)$ is isometric (up to a constant multiple) to $(M_1,g_1)$.
Finally, if neither $u_1$ nor $u_2$ is constant, then, $(M_1,g_1)$ and $(M_2,g_2)$ are both round $n$-spheres, so $M_0$ is as well.