In the absence of more clever ideas, I like using coordinates. Without loss of generality you can pick three points as a projective basis on the line. This will reduce the number of variables somewhat. Pretty much at random I did pick
$$A'=\begin{pmatrix}1\\0\end{pmatrix}\qquad
B'=\begin{pmatrix}0\\1\end{pmatrix}\qquad
U=\begin{pmatrix}1\\1\end{pmatrix}$$
Pick one coordinate for one of the other points:
$$P=\begin{pmatrix}p\\1\end{pmatrix}$$
The rest can be derived from these. A harmonic throw has cross-ratio $-1$. Writing $[\cdot,\cdot]$ for a $2\times2$ determinant the harmonic throw $APUA'$ translates into
\begin{align*}
-1=(A,P;U,A')&=\frac{[A,U][P,A']}{[A,A'][P,U]}\\
-[A,A'][P,U]&=[A,U][P,A']\\
0&=[A,U][P,A']+[A,A'][P,U]\\
0&=[A,[P,A']U+[P,U]A']\\
A&\sim[P,A']U+[P,U]A'\\
A&\sim\begin{vmatrix}p&1\\1&0\end{vmatrix}\begin{pmatrix}1\\1\end{pmatrix}+
\begin{vmatrix}p&1\\1&1\end{vmatrix}\begin{pmatrix}1\\0\end{pmatrix}\\
A&=\begin{pmatrix}2-p\\1\end{pmatrix}
\end{align*}
Of course you could have taken any other representative, but this one is particularly convenient. In the same fashion from $(B,P;U,B')=-1$ you get
$$[P,B']U+[P,U]B'=
\begin{vmatrix}p&0\\1&1\end{vmatrix}\begin{pmatrix}1\\1\end{pmatrix}+
\begin{vmatrix}p&1\\1&1\end{vmatrix}\begin{pmatrix}0\\1\end{pmatrix}=
\begin{pmatrix}p\\2p-1\end{pmatrix}=B$$
And from $(A',B';P,Q)=-1$ you get
$$[P,A']B'+[P,B']A'=
\begin{vmatrix}p&1\\1&0\end{vmatrix}\begin{pmatrix}0\\1\end{pmatrix}+
\begin{vmatrix}p&0\\1&1\end{vmatrix}\begin{pmatrix}1\\0\end{pmatrix}\sim
\begin{pmatrix}-p\\1\end{pmatrix}=Q$$
Finally from $(A,B;U,O)=-1$ you get (simplified by cancelling a factor of $2p-2$)
$$[U,A]B+[U,B]A=
\begin{vmatrix}1&2-p\\1&1\end{vmatrix}\begin{pmatrix}p\\2p-1\end{pmatrix}+
\begin{vmatrix}1&p\\1&2p-1\end{vmatrix}\begin{pmatrix}2-p\\1\end{pmatrix}\sim
\begin{pmatrix}1\\p\end{pmatrix}=O$$
Given these points, you can formulate the transformation matrices explicitly using this computation.
$$\Lambda=\begin{pmatrix}1-2p&p\\p&p^2-2p\end{pmatrix}\qquad
\Theta=\begin{pmatrix}2p-1&-p\\1&p-2\end{pmatrix}$$
Now it's a matter of doing the multiplications to confirm that
$$[\Lambda P,\Theta O]=0$$
i.e. the determinant is zero as the vectors are scalar multiples of one another. For the inverse operations, you can use the adjugate matrix, as this matches the inverse up to some factor.
$$\operatorname{adj}\Lambda=\begin{pmatrix}p^2-2p&-p\\-p&1-2p\end{pmatrix}\qquad
\operatorname{adj}\Theta=\begin{pmatrix}p-2&p\\-1&2p-1\end{pmatrix}$$
With these you can check the second half of the claim in the same fashion.
