This is a coordinate-based proof making heavy use of computer algebra. It'll only answer the first of your questions.
Parametrization
Your problem is invariant under Möbius transformations, since it only deals with incidences between circles and points. A Möbius transformation is uniquely defined by three points and their images. So without loss of generality, I choose $A=(-1,0),B=(0,1),D=(1,0)$ so that the circle $ABCDEF$ will be the unit circle. Using the tangent half-angle substitution, you can choose $C=(1-c^2,2c)/(1+c^2)$ and likewise for $E$ and $F$ defined by one real parameter $e$ resp. $f$ each. The point $A_1=(x,y)$ can be chosen arbitrarily, and that choice fixes the circle $C_{AD}$. As $C_{BE}$ is any circle through $B$ and $E$, it must be a linear combination of the unit circle and the line $BE$. So I write $C_{BE}=BE+m\cdot U$ where $U$ denotes the unit circle $ABCDEF$ and $m$ is any real parameter. Likewise $C_{CF}=CF+n\cdot U$. After this, everything is fixed. So the whole setup has $7$ real degrees of freedom, modeled by the variables $c,e,f,x,y,m,n$ in my parametrization.
Möbius geometry
For most of my computation, I'll be using $\mathbb P(\mathbb R^{1,3})$ to represent circles, lines and points in a uniform way. $\mathbb R^{1,3}$ means a four-dimensional vector space with an inner product that has signature $(+,-,-,-)$, and $\mathbb P$ indicates that it's a projective space using homogeneous coordinates, so multiples of a vector will represent the same object. I'll be describing a circle with center $x,y$ and radius $r$ as a vector
$$\begin{pmatrix}x\\y\\x^2+y^2-r^2\\1\end{pmatrix}$$
or any multiple thereof. To match this format, I'll define the inner product as $\langle a,b\rangle:=a^T\cdot M\cdot b$ with
$$M=\begin{pmatrix}-2&0&0&0\\0&-2&0&0\\0&0&0&1\\0&0&1&0\end{pmatrix}$$
Others often choose a different basis, so that the matrix of the inner product only has $\pm 1$ on the diagonal. But I prefer separating the vector elements the way I did. Circles $a$ and $b$ intersect orthogonally iff $\langle a,b\rangle=0$, but we won't need this general property here. A vector $p$ is a point (i.e. a circle with zero radius) iff $\langle p,p\rangle=0$. A point $p$ lies on a circle $c$ iff $\langle p,c\rangle=0$. In all of this, lines are just special cases of circles. A line $aX+bY=c$ would correspond to the vector $(a,b,2c,0)^T$.
With the parametrization above you'd have circles like these:
$$
A=\begin{pmatrix}-1\\0\\1\\1\end{pmatrix}\quad
B=\begin{pmatrix}0\\1\\1\\1\end{pmatrix}\quad
C=\begin{pmatrix}1-c^2\\2c\\1+c^2\\1+c^2\end{pmatrix}\quad
A_1=\begin{pmatrix}x\\y\\x^2+y^2\\1\end{pmatrix}\quad
U=\begin{pmatrix}0\\0\\-1\\1\end{pmatrix}
$$
Note how I've moved the denominator from the original definition of $C$ into a different choice of representatives. That way I can avoid divisions, which is an important goal not only in this step but also in later steps, since dealing with polynomials is often easier than dealing with rational functions.
Circle through three points
Using this formalism, the task of finding a circle through three given points $a,b,c$ is equivalent to finding a vector $v$ which satisfies $a^T\cdot M\cdot v=b^T\cdot M\cdot v=c^T\cdot M\cdot v=0$. This is an underspecified linear system of equations, which in general will have a one-dimensional space of solutions representing the circle in question. One way to obtain an element of that space is by writing the vectors $Ma$ through $Mc$ into a $3\times4$ matrix $K$. Then compute the minors $K_i$ as the determinant of the submatrix you get when you drop column $i$ from the matrix. $(-K_1,K_2,-K_3,K_4)^T$ will be one element of the kernel of that matrix. So for example $C_{AD}$ is the circle through $A,D,A_1$:
$$
K=\begin{pmatrix}
2 & 0 & 1 & 1 \\
-2 & 0 & 1 & 1 \\
-2x & -2y & 1 & x^2+y^2
\end{pmatrix}
\qquad
C_{AD}
=\begin{pmatrix} 0 \\ 4(x^2+y^2-1) \\ -8y \\ 8y \end{pmatrix}
\sim\begin{pmatrix} 0 \\ x^2+y^2-1 \\ -2y \\ 2y \end{pmatrix}
$$
This same algorithm can also be used to find the line joining two points. That's because lines are just special circles passing through the point at infinity, which has coordinates $(0,0,1,0)^T$. So the circle passing through two finite points and the line at infinity will represent the line joining the finite points. One can use this to describe e.g. $BE$, the line joining $B$ and $E$, which I've used to represent $C_{BE}=BE+m\cdot U$.
Other point of intersection
The second essential operation I've used is the following: Given tow circles $a$ and $b$ and a point $p$ incident with both of them, find the other point of intersection.
Starting with the two circle vectors $a$ and $b$, the first step is finding the radical line of these two circles, i.e. the line passing through both their points of intersection. All linear combinations $\lambda a+\mu b$ must pass through these points of intersection because $\langle\lambda a+\mu b,p\rangle=\lambda\langle a,p\rangle+\mu\langle b,p\rangle$ so if the latter two inner products are zero, the linear combination can be nothing else. The radical line is the uniquely defined linear combination which also passes through the point at infinity, i.e. satisfies $\langle\lambda a+\mu b,p_\infty\rangle=0$. In order to elegantly describe that, choose $\lambda:=\langle b,p_\infty\rangle$ and $\mu:=-\langle a,p_\infty\rangle$. That trick is sometimes known as Plücker's µ. So $l=\langle b,p_\infty\rangle\cdot a-\langle a,p_\infty\rangle\cdot b$ is the radical line. As a line, it is of the form $(r_1,r_2,r_3,0)^T$ representing the equation $r_1X+r_2Y=r_3/2$.
At this point I decided to drop the Möbius formalism for a moment, and instead move to the formalism used for conics in real projective geometry. That's because in Möbius geometry, there only exists a single point at infinity, whereas in the real projective plane there is a whole line at infinity, and the radical line would intersect the line at infinity at a point $f=(r_2,-r_1,0)^T$.
In order to represent a circle with center $(x,y)$ and radius $r$ as a conic, observe that it corresponds to the set of all points $(X,Y)$ satisfying
$$0=(X-x)^2+(Y-y)^2-r^2= %X^2-2xX+x^2+Y^2-2yY+y^2-r^2=
(X,Y,1)\cdot
\begin{pmatrix}1&0&-x\\0&1&-y\\-x&-y&x^2+y^2-r^2\end{pmatrix}
\cdot\begin{pmatrix}X\\Y\\1\end{pmatrix}$$
That symmetric matrix is used to represent the circle as a conic. It is homogeneous, too, since multiples of the equation describe the same set of points. Comparing it with the circle vectors introduced above, it is easy to see that the vector $a$ corresponds to a matrix
$$a'=\begin{pmatrix}a_4&0&-a_1\\0&a_4&-a_2\\-a_1&-a_2&a_3\end{pmatrix}$$
Any point on a line can be represented as a linear combination of two distinct points spanning that line. So the line through $p'=(p_1,p_2,p_4)^T$ and $f$ is the set of points $q'=\lambda p'+\mu f$. In order for a point to lie on that line and on the circle $a$ (i.e. the conic $a'$), that linear combination has to satisfy
$$0=q'^T\cdot a'\cdot q=
\lambda^2(p'^T\cdot a'\cdot p')+
2\lambda\mu(p'^T\cdot a'\cdot f)+
\mu^2(f^T\cdot a'\cdot f)$$
which at first looks quadratic. But knowing $p'^T\cdot a'\cdot p'=0$ (since $p$ was assumed to lie on both circles), one can cancel one $\mu$, solve the resulting homogeneous equation using Plücker's µ again, and write
$$q'=(f^T\cdot a'\cdot f)\cdot p' - 2(p'^T\cdot a'\cdot f)\cdot f$$
This corresponds to a radius zero circle $({q'}_1{q'}_3,{q'}_2{q'}_3,{q'}_1^2+{q'}_2^2,{q'}_3^2)^T$ in order to avoid division by ${q'}_3$.
Putting it all together
So now one has the two essential algorithms, namely circle through three points and other intersection, expressed as algorithms who only perform ring operations on the individual coordinates, no divisions. So starting from the parametrization, one can perform the whole computation using vectors of polynomials at each step.
It helps a lot to cancel common factors after each such computation. Those common factors would correspond to situations where some intermediate result becomes undefined. Canceling them corresponds to lifting a liftable singularity. And it helps to keep the size of the polynomials down, to make computations feasible.
As the result of such a chain of computations, one would end up with coordinates for all the points and circles involved. Not particularly simple polynomials; I've seen coordinates of degree up to $8$ for some of the points, although the circles (e.g. through $B,C,B_1,C_1$) are simpler with a maximal degree of $4$. I won't paste those vectors here since they take up a lot of room and add little value.
In order to check the cocircularity of $A_1,A,F,F_1$, one could simply compute the determinant of the four corresponding Möbius vectors. It is zero, indicating the existence of a circle (not the null vector) incident with all four of them. For the cocircularity of $A_1,B_1,C_1,D_1,E_1,F_1$ one could either do multiple determinants, or compute the circle through $A_1,B_1,C_1$ as described above, and then check that it's incident with the other three points by computing the inner product. By the way, that circle has coordinates of degree up to $6$.
I've done all of the above using the Sage code I posted here, and could verify your claim that way. It would be possible, though tedious, to print all the vectors and verify all the incidence conditions manually. So there is nothing here which means that this proof needs a computer, but doing it without one takes considerable work and at least I'm too lazy for that. Looking forward to other approaches, though.
