Write $C=-A$, $D=-B$ and $F=-E$ (so their common midpoint is $0$) and define \begin{align*} P&=+A+E-B\\ Q&=-A+E+B\\ R&=-A-E-B\\ S&=+A-E+B \end{align*}It is easy to check that $A$, $B$, $C$, $D$, $E$ und $F$ are the midpoints of $PS$, $QS$, $QR$, $PR$, $PQ$ and $RS$. Then the $4$ circles in the problem are the nine-point-circles of $PQS$, $QRS$, $PQR$ and $PRS$, which concur at the Poncelet-Point of $PQRS$ (or are all the same if $P$, $Q$, $R$ and $S$ form an orthocentric system). sketch[asy][asy] /* Geogebra to Asymptote conversion, documentation at artofproblemsolving.com/Wiki go to User:Azjps/geogebra */ import graph; size(12cm); real labelscalefactor = 0.5; /* changes label-to-point distance */ pen dps = linewidth(0.7) + fontsize(10); defaultpen(dps); /* default pen style */ pen dotstyle = black; /* point style */ real xmin = -10., xmax = 10., ymin = -10., ymax = 8.; /* image dimensions */ pair O = (0.,0.), A = (3.006452682926829,-4.600341885837027), C = (-3.006452682926829,4.600341885837027), B = (0.5496214634146326,-2.029239210723372), D = (-0.5496214634146326,2.029239210723372), F = (-3.0635882926829265,-2.2844997219673155), P = (5.520419512195122,-0.28660295314633943), Q = (0.6067570731707299,4.85560239708097), R = (-6.619662439024388,4.345081374593084); /* draw figures */ draw(D--B, linewidth(1.2)); draw(F--(3.0635882926829265,2.2844997219673155), linewidth(1.2)); draw(C--A, linewidth(1.2)); draw((0.49248585365853526,-8.914080818527715)--P, linewidth(1.2)); draw(Q--(0.49248585365853526,-8.914080818527715), linewidth(1.2)); draw(Q--P, linewidth(1.2)); draw(R--Q, linewidth(1.2)); draw(R--(0.49248585365853526,-8.914080818527715), linewidth(1.2)); draw(R--P, linewidth(1.2)); /* dots and labels */ dot(O,dotstyle); label("$O$", (0.06396878048780306,0.48472784938775715), NE * labelscalefactor); dot(A,dotstyle); label("$A$", (3.4921053658536585,-4.686045308340816), NE * labelscalefactor); dot(C,dotstyle); label("$C$", (-2.878515121951223,4.884170204582234), NE * labelscalefactor); dot(B,dotstyle); label("$B$", (0.778163902439023,-2.429188515741052), NE * labelscalefactor); dot(D,dotstyle); label("$D$", (-0.4216839024390265,2.3130675294685785), NE * labelscalefactor); dot((3.0635882926829265,2.2844997219673155),dotstyle); label("$E$", (3.1778595121951216,2.570177796979944), NE * labelscalefactor); dot(F,dotstyle); label("$F$", (-3.021354146341467,-1.743561135710744), NE * labelscalefactor); dot(P,dotstyle); label("$P$", (5.634690731707318,-9.248781337110491e-4), NE * labelscalefactor); dot((0.49248585365853526,-8.914080818527715),dotstyle); label("$S$", (0.9781385365853645,-8.999784241031504), NE * labelscalefactor); dot(Q,dotstyle); label("$Q$", (0.7210282926829253,5.141280472093599), NE * labelscalefactor); dot(R,dotstyle); label("$R$", (-6.50662634146342,4.627059937070868), NE * labelscalefactor); clip((xmin,ymin)--(xmin,ymax)--(xmax,ymax)--(xmax,ymin)--cycle); /* end of picture */ [/asy][/asy]

Again define $c=-a$, $d=-b$ and $f=-e$ and let $z$ be the second intersection of $(AEB)$ and $(BFC)$, assuming they're not already the same circle (in which case by PoP on $0$ we have $|a|=|e|$ and thus $0$ is the center of this circle. Thus $D$ must also lie on the Circle and we're done). Then we know that $\frac{z-b}{z-a}\cdot\frac{e-a}{e-b}\in\mathbb R$ and $\frac{z-b}{z-f}\cdot\frac{c-f}{c-b}\in\mathbb R$. Dividing, we get \[\mathbb R\ni\frac{z-b}{z-b}\cdot\frac{z-f}{z-a}\cdot\frac{e-a}{c-f}\cdot\frac{c-b}{e-b}=\frac{z-f}{z-a}\cdot\frac{c-b}{e-b}=\frac{z-f}{z-a}\cdot\frac{d-a}{d-f}\]by repeatedly using the definitions above. But that just means that $ZFDA$ is cyclic as well, as desired. But what happens if $z=b$, i.e. circles $(AEB)$ and $(BFC)$ are tangent in $B$? We could, for example, reformulate the entire proof in terms of the angles the lines form and get a very quick angle chase where line $BZ$ is just the common tangent, and everything works out again (in fact, the "complex-bash"-proof above is really just a reformulation of an angle chase proof)