[proofplan]
We prove the identity by writing both matrices in coordinates and comparing the lower row of the product matrix with the product of the two scalar factors. The only algebraic point is that the denominator appearing in $\gamma_2 z$ is exactly $j(\gamma_2,z)$, so it cancels when the two factors are multiplied. Since the automorphy factor depends only on the lower row of the matrix, equality of the resulting linear expressions in $z$ proves the cocycle identity.
[/proofplan]
[step:Write the two matrices and compute the lower row of their product]
Let
\begin{align*}
\gamma_1 =
\begin{pmatrix}
a_1 & b_1 \\
c_1 & d_1
\end{pmatrix},
\qquad
\gamma_2 =
\begin{pmatrix}
a_2 & b_2 \\
c_2 & d_2
\end{pmatrix}
\end{align*}
be elements of $SL_2(\mathbb{Z})$. Matrix multiplication gives
\begin{align*}
\gamma_1\gamma_2
&=
\begin{pmatrix}
a_1a_2+b_1c_2 & a_1b_2+b_1d_2 \\
c_1a_2+d_1c_2 & c_1b_2+d_1d_2
\end{pmatrix}.
\end{align*}
By the definition of $j$, which uses the lower row of the matrix, we therefore have
\begin{align*}
j(\gamma_1\gamma_2,z)
&=
(c_1a_2+d_1c_2)z+(c_1b_2+d_1d_2).
\end{align*}
[/step]
[step:Expand the product of the two automorphy factors]
First $c_2z+d_2\neq 0$. Indeed, if $c_2=0$, then the determinant condition $a_2d_2-b_2c_2=1$ gives $a_2d_2=1$, so $d_2\neq 0$. If $c_2\neq 0$ and $c_2z+d_2=0$, then $z=-d_2/c_2\in\mathbb{R}$, contradicting $z\in\mathbb{H}$. Thus the fractional linear action is defined at $z$, and
\begin{align*}
\gamma_2 z = \frac{a_2z+b_2}{c_2z+d_2}.
\end{align*}
Moreover $\gamma_2z\in\mathbb{H}$ because $\gamma_2\in SL_2(\mathbb{R})$ and
\begin{align*}
\operatorname{Im}(\gamma_2z)
&= \operatorname{Im}\left(\frac{a_2z+b_2}{c_2z+d_2}\right)\\
&= \frac{\operatorname{Im}(z)}{|c_2z+d_2|^2} > 0,
\end{align*}
where the displayed identity follows by multiplying numerator and denominator by the complex conjugate of $c_2z+d_2$ and using $a_2d_2-b_2c_2=1$. Hence $\gamma_2 z$ is in the domain of $j(\gamma_1,\cdot)$ and $z$ is in the domain of $j(\gamma_2,\cdot)$. The two factors are
\begin{align*}
j(\gamma_1,\gamma_2 z)
&= c_1(\gamma_2 z)+d_1
= c_1\frac{a_2z+b_2}{c_2z+d_2}+d_1,\\
j(\gamma_2,z)
&= c_2z+d_2.
\end{align*}
Multiplying these expressions and collecting terms gives
\begin{align*}
j(\gamma_1,\gamma_2 z)j(\gamma_2,z)
&=
\left(c_1\frac{a_2z+b_2}{c_2z+d_2}+d_1\right)(c_2z+d_2)\\
&=
c_1(a_2z+b_2)+d_1(c_2z+d_2)\\
&=
(c_1a_2+d_1c_2)z+(c_1b_2+d_1d_2).
\end{align*}
[/step]
[step:Compare the two expressions]
The expression obtained for $j(\gamma_1\gamma_2,z)$ is
\begin{align*}
(c_1a_2+d_1c_2)z+(c_1b_2+d_1d_2),
\end{align*}
and the expression obtained for $j(\gamma_1,\gamma_2 z)j(\gamma_2,z)$ is the same. Hence
\begin{align*}
j(\gamma_1\gamma_2,z)=j(\gamma_1,\gamma_2 z)j(\gamma_2,z).
\end{align*}
This proves the cocycle identity for all $\gamma_1,\gamma_2 \in SL_2(\mathbb{Z})$ and all $z \in \mathbb{H}$.
[/step]