Find a fractional linear transformation $f$ that maps $\{z:\ Im(z)>0\}$ onto $\{z: |z|<1\}$ and such that $f(i)=0$ and $Re (f'(i))=0.$
I consider $x$-axis maps to the unit circle. First, let $f(z)=(az+b)/(cz+d)$. Since $f(i)=0$, I can obtain that $b=-ai$. Then, I tried to choose points from $x-$axis to unit circle to determine the function, let $f(0)=1,\ f(1)=i,\ f(-1)=-i.$ By some computation, I get $$ f(z)=\frac{az-ai}{-az-ai}=\frac{-z+i}{z+i}. $$ (The singularity is $-i$, and the $x-$axis doesn't pass through the point $-i$, so it guarantees that the line with maps to a circle.) I check with the derivative: $$ f'(z)=\frac{-(z+i)-(-z+i)}{(z+i)^2}=\frac{-2i}{(z+i)^2}. $$ Then, it is clear that $Re (f'(i))=0$.
For me, the function I found seems works to the problem. But, I wonder is there any standard way to deal with this problem? For my feeling, it is just lucky that what I choose ($f(0)=1,\ f(1)=i,\ f(-1)=-i$) satisfy the condition $Re(f'(i))=0$.
The equation you stumbled upon is actually quite well known (for people that do this kind of thing) and is the Cayley Transform.
Really, the way you did it is the way you want to think of it. Start with the general Mobius Map and deduce what the values need to be. But I'm guessing you want a conceptual/heuristic motivation for how someone would think to do this?
To this end think about what we're trying to do. We want to map a half-plane onto a unit circle. (Note: we can think about this in two ways, either circle -> half plane or half plane -> circle).
Recall that in complex analysis circles and lines are essentially the same thing (as explained in this stack exchange response). So, we need to figure out what circle and what line make sense to match up.
But to do this all we really need to think about is the boundary. For the circle we want, the boundary is $|z| = 1$ (a circle). For the half plane we want, the boundary is $Im(z) = 0$, so we want to map one of these onto the other (depending on which way we want to go). (Note: This is a heuristic approach. Technically those aren't the entire boundary, but they are the obvious pieces to think about as they are the 'restrictions' that we imposed).
This should rather quickly lead you to the map sending 3 points on the imaginary axis to 3 points on the unit circle. Which three don't make a huge difference initially, as long as you make sure you orient it correctly. Once you get some map from the axis to the circle, it's a matter of tweaking the coefficients to find the one that conforms to your given criteria (i.e. the specific value of the derivative and the specific point of the transform), which is the rabbit hole you then went down.
Hope that helps. If I completely misinterpreted your question, feel free to let me know and I'll try to adjust.