I happen to have some notes on this question. What follows here is the usual approach, it's just multivariate calculus paired with the Cauchy Riemann equations. I have an idea for an easier way, I'll post it as a second answer in a bit if it works.
If we use polar coordinates to rewrite $f$ as follows:
$$ f(x(r,\theta),y(r,\theta)) = u(x(r,\theta),y(r,\theta))+iv(x(r,\theta),y(r,\theta)) $$
we use shorthands $F(r,\theta)=f(x(r,\theta),y(r,\theta))$ and $U(r,\theta )=u(x(r,\theta),y(r,\theta))$ and $V(r,\theta )=v(x(r,\theta),y(r,\theta))$. We derive the CR-equations in polar coordinates via the chain rule from multivariate calculus,
$$ U_r = x_ru_x + y_ru_y = \cos(\theta)u_x + \sin(\theta)u_y \ \
\text{and} \ \ U_{\theta} = x_{\theta}u_x + y_{\theta}u_y = -r\sin(\theta)u_x + r\cos(\theta)u_y $$
Likewise,
$$ V_r = x_rv_x + y_rv_y = \cos(\theta)v_x + \sin(\theta)v_y \ \
\text{and} \ \ V_{\theta} = x_{\theta}v_x + y_{\theta}v_y = -r\sin(\theta)v_x + r\cos(\theta)v_y $$
We can write these in matrix notation as follows:
$$ \left[ \begin{array}{l} U_r \\ U_{\theta} \end{array} \right] = \left[ \begin{array}{ll} \cos(\theta) & \sin(\theta) \\ -r\sin(\theta) & r\cos(\theta) \end{array} \right]\left[ \begin{array}{l} u_x \\ u_y \end{array} \right] \ \ \text{and} \ \
\left[ \begin{array}{l} V_r \\ V_{\theta} \end{array} \right] = \left[ \begin{array}{ll} \cos(\theta) & \sin(\theta) \\ -r\sin(\theta) & r\cos(\theta) \end{array} \right]\left[ \begin{array}{l} v_x \\ v_y \end{array} \right] $$
Multiply these by the inverse matrix: $\left[ \begin{array}{ll} \cos(\theta) & \sin(\theta) \\ -r\sin(\theta) & r\cos(\theta) \end{array} \right]^{-1} = \frac{1}{r}\left[ \begin{array}{ll} r\cos(\theta) & -\sin(\theta) \\ r\sin(\theta) & \cos(\theta) \end{array} \right]$ to find
$$ \left[ \begin{array}{l} u_x \\ u_y \end{array} \right] = \frac{1}{r}\left[ \begin{array}{ll} r\cos(\theta) & -\sin(\theta) \\ r\sin(\theta) & \cos(\theta) \end{array} \right]\left[ \begin{array}{l} U_r \\ U_{\theta} \end{array} \right] = \left[ \begin{array}{l} \cos(\theta)U_r - \tfrac{1}{r}\sin(\theta)U_{\theta} \\
\sin(\theta)U_r + \tfrac{1}{r}\cos(\theta)U_{\theta} \end{array} \right] $$
A similar calculation holds for $V$. To summarize:
$$ u_x = \cos(\theta)U_r - \tfrac{1}{r}\sin(\theta)U_{\theta} \ \ \ \ v_x = \cos(\theta)V_r - \tfrac{1}{r}\sin(\theta)V_{\theta} $$
$$ u_y =\sin(\theta)U_r + \tfrac{1}{r}\cos(\theta)U_{\theta} \ \ \ \ v_y =\sin(\theta)V_r + \tfrac{1}{r}\cos(\theta)V_{\theta} $$
The CR-equation $u_x=v_y$ yields:
$$ (A.) \ \ \cos(\theta)U_r - \tfrac{1}{r}\sin(\theta)U_{\theta} = \sin(\theta)V_r + \tfrac{1}{r}\cos(\theta)V_{\theta} $$
Likewise the CR-equation $u_y=-v_x$ yields:
$$ (B.) \ \ \sin(\theta)U_r + \tfrac{1}{r}\cos(\theta)U_{\theta} = -\cos(\theta)V_r + \tfrac{1}{r}\sin(\theta)V_{\theta}$$
Multiply (A.) by $r\sin(\theta)$ and $(B.)$ by $r\cos(\theta)$ and subtract (A.) from (B.):
$$ \boxed{U_{\theta} = -rV_r} $$
Likewise multiply (A.) by $r\cos(\theta)$ and $(B.)$ by $r\sin(\theta)$ and add (A.) and (B.):
$$ \boxed{rU_r = V_{\theta}} $$
Finally, recall that $z = re^{i\theta}=r(\cos(\theta)+i\sin(\theta))$ hence
\begin{align} \notag
f'(z) &= u_x+iv_x \\ \notag
&= (\cos(\theta)U_r - \tfrac{1}{r}\sin(\theta)U_{\theta})+i(\cos(\theta)V_r - \tfrac{1}{r}\sin(\theta)V_{\theta}) \\ \notag
&= (\cos(\theta)U_r + \sin(\theta)V_{r})+i(\cos(\theta)V_r - \sin(\theta)U_{r}) \\ \notag &= (\cos(\theta)- i\sin(\theta))U_r + i(\cos(\theta)-i\sin(\theta))V_r \\ \notag
&= e^{-i\theta}( U_r+iV_r) \notag
\end{align}
Here's the Taylor expansion approach. First, we need a linear algebra fact:
Any $\mathbb R$-linear map $T:\mathbb C\to\mathbb C$ can be written as $Tz=\alpha z+\beta\bar z$ for some (unique) $\alpha,\beta\in \mathbb C$.
Uniqueness should be clear, and existence follows by first writing $T(x+iy)=\gamma x+\delta y$ with complex $\gamma,\delta$ and then replacing $x=(z+\bar z)/2$, $y=(z-\bar z)/(2i)$.
This linear algebra fact is worth remembering, as it simplifies various computations in complex analysis.
Back to the problem. The first-order real Taylor expansion of $f$ at a point $z\in\mathbb C$ takes the form
$$
f(z+h)=f(z)+\alpha h+\beta \bar h+o(|h|)
$$
and we want to show that $\beta=0$ here. The equations you are given say that $r f_r + i f_\theta = 0$.
Write $h=r e^{i\theta}$, and take the derivatives of $\alpha h+\beta \bar h$:
$$
\begin{align}
f_r &= \alpha e^{i\theta} + \beta e^{-i\theta} \\
f_\theta & = \alpha ir e^{i\theta} - \beta i r e^{-i\theta} \\
rf_r+f_\theta & = 2\beta r e^{-i\theta}
\end{align}
$$
Thus $\beta=0$ as desired.
Best Answer
Take the Cauchy-Riemann equations in polar form, that is $$u_r = \frac{1}{r}v_{\theta},\quad \frac{1}{r} u_{\theta} = -v_r .$$ Now, take the partial derivative of the first equation with respect to $r$ and the partial derivative of the second equation with respect to $\theta$, then $$ u_{rr} = \frac{\partial}{\partial r} (\frac{1}{r}v_{\theta}) = \frac{v_{\theta\ r}}{r} - \frac{v_{\theta}}{r^2},\\ u_{\theta \theta} = \frac{\partial}{\partial \theta}(-r v_r)= - r v_{r \theta} .$$ Here, I just took the partial derivative for the first one with respect to $r$ and the partial derivative with respect to $\theta$ for the second one, but you can work out the other two partial derivatives. Now, assuming continuity of partial derivatives, one can say that $v_{r \theta} = v_{\theta r},$ so we can substitute the last equation on the one before, thus $$u_{rr} = -\frac{u_{\theta \theta}}{r^2} - \frac{v_{\theta}}{r^2},$$ but we can again make use of the Cauchy-Riemann equations and substitue $v_{\theta} = r u_r,$ and using this result in the last equation yields $$ r^2 u_{rr} = -u_{\theta \theta} - r u_r,\\r^2 u_{rr} + r u_r + u_{\theta \theta} = 0.$$