At some point you have $b=-c(x^2+y^2)$ and $a=2cx+d$. If $c=0$, then $b=0$, in which case the element $\begin{pmatrix} a&b\\c&d\end{pmatrix}$ of the projective special linear group $\operatorname{PSL}_2(\Bbb Z)$ of this form is the identity element (recalling that $\operatorname{PSL}_2(\Bbb Z)=\operatorname{SL}_2(\Bbb{Z})/\{\pm I\}$), so we are not interested in this case. So, $c$ is assumed to be non-zero, and can without loss of generality assumed to be positive.
That is, $x=\frac{a-d}{2c}$, so that $$-\frac{b}{c}=x^2+y^2=\frac{(a-d)^2}{4c^2}+y^2.$$
Since $y>0$, we have
$$y=\frac{\sqrt{4-(a+d)^2}}{2c},$$ which means $d=-a$ or $d=-a\pm1$.
In the case $d=-a$, we must have $x=\frac{a}{c}$ and $y=\frac{1}{c}$.
For a given point $z=x+yi$ of this form that is fixed by a non-trivial $\gamma\in\operatorname{PSL}_2(\Bbb{Z})$, we can assume that $$\gamma=\begin{pmatrix}a&b\\c&-a\end{pmatrix}\wedge bc=-(a^2+1).\tag{1}$$ That is, $\gamma^2$ is the identity of $\operatorname{PSL}_2(\Bbb{Z})$. The points in $\Bbb H$ with a stabilizer of order $2$ are of the form $z=x+yi$, where
$$x=\frac{a}{c}\wedge y=\frac{1}{2c}$$
for some integer $a$ and for some integer $c>0$, such that $a^2+1$ is divisible by $c$ (so that there exists $b\in \Bbb Z$ such that $bc=-(a^2+1)$). The number of orbits for a given $a$ is $\sigma_0(a^2+1)$, where $\sigma_0$ is the divisor counting function.
For the case $d=-a+1$, we see that $x=\frac{2a-1}{2c}$ and $y=\frac{\sqrt{3}}{2c}$. Then, $z=x+yi$ is fixed by
$$\gamma=\begin{pmatrix}a&b\\c&-a+1\end{pmatrix}\wedge bc=-(a^2-a+1).\tag{2}$$
Note that $\gamma^3$ is the identity of $\operatorname{PSL}_2(\Bbb{Z})$. Thus, for a given $a$, there are corresponding $\sigma_0(a^2-a+1)$ points $z\in \Bbb H$.
For the case $d=-a-1$, we see that $x=\frac{2a+1}{2c}$ and $y=\frac{\sqrt{3}}{2c}$. Then, $z=x+yi$ is fixed by
$$\gamma=\begin{pmatrix}a&b\\c&-a-1\end{pmatrix}\wedge bc=-(a^2+a+1).\tag{3}$$
Note that $\gamma^3$ is the identity of $\operatorname{PSL}_2(\Bbb{Z})$. Thus, for a given $a$, there are corresponding $\sigma_0(a^2+a+1)$ points $z\in \Bbb H$. (It can be seen that this case is identical to the previous case via the transformation $a\mapsto a-1$.)
For example, $x=\frac{3}{5}$ and $y=\frac{1}{5}$ fit the bill (with $a=3$, $b=-2$, $c=5$, and $d=-a=-3$) for a point with stabilizer of order $2$, with $$\gamma=\begin{pmatrix}3&-2\\5&-3\end{pmatrix}.$$
For the same $a=3$, there are three more points with stablizers of order $2$, i.e., with $c=1$, $c=2$, and $c=10$. That is, for $a=3$, we have in total four points with non-trivial stabilizers $\gamma$ of order $2$: $z=3+i$ with $\gamma=\begin{pmatrix}3&-10\\1&-3\end{pmatrix}$, $z=\frac{3}{2}+\frac{i}{2}$ with $\gamma=\begin{pmatrix}3&-5\\2&-3\end{pmatrix}$, $z=\frac{3}{5}+\frac{i}{5}$ with $\gamma=\begin{pmatrix}3&-2\\5&-3\end{pmatrix}$, and $z=\frac{3}{10}+\frac{i}{10}$ with $\gamma=\begin{pmatrix}3&-1\\10&-3\end{pmatrix}$.
For the same $a=3$, there are also four points with stabilizers of order $3$. That is, for $a=3$, we have four points with non-trivial stabilizers $\gamma$ of order $3$: $z=\frac52+\frac{\sqrt{3}i}{2}$ with $\gamma=\begin{pmatrix}3&-7\\1&-2\end{pmatrix}$, $z=\frac{5}{14}+\frac{\sqrt{13}i}{14}$ with $\gamma=\begin{pmatrix}3&-1\\7&-2\end{pmatrix}$, $z=\frac{7}{2}+\frac{\sqrt{3}i}{2}$ with $\gamma=\begin{pmatrix}3&-13\\1&-4\end{pmatrix}$, and $z=\frac{7}{26}+\frac{\sqrt{3}i}{26}$ with $\gamma=\begin{pmatrix}3&-1\\13&-4\end{pmatrix}$.
In conclusion, there are three kinds of points $z\in\Bbb H$---those with trivial stabilizers, those with stabilizers of order $2$, and those with stabilizers of order $3$. The stablizers in non-trivial cases are generated by $\gamma$ given in (1) and (2).
