Deriving the equations of motion of a string from Polyakov action The Polyakov action  for a string in spacetime $(M, g)$ is
$$\mathcal{L} = \frac{T}{2}\sqrt{-\det(h)}h^{\eta\lambda}(\sigma)g_{\mu\nu}(x(\sigma))\frac{\partial x^{\mu}}{\partial \sigma^{\eta}}\frac{\partial x^{\nu}}{\partial \sigma^{\lambda}}$$
And Euler-Lagrange equations for a classical field reads
$$\frac{\partial}{\partial \sigma^{\rho}}\frac{\partial \mathcal{L}}{\partial(\frac{x^{\alpha}}{\partial \sigma^{\rho}})}-\frac{\partial \mathcal{L}}{\partial x^{\alpha}} = 0$$
Using the fact that $\displaystyle \frac{\partial \mathcal{L}}{\partial x^{\alpha}} = \frac{T}{2}\sqrt{-\det(h)}h^{\eta\lambda}(\sigma)\frac{\partial g_{\mu\nu}}{\partial x^{\alpha}}\frac{\partial x^{\mu}}{\partial \sigma^{\eta}}\frac{\partial x^{\nu}}{\partial \sigma^{\lambda}}$, we conclude that
$$\frac{\partial}{\partial \sigma^{\rho}}\biggr(\frac{T}{2}\sqrt{-\det(h)}h^{\eta\rho}(\sigma)g_{\mu\alpha}(x(\sigma))\frac{\partial x^{\mu}}{\partial \sigma^{\eta}}\biggr) - \frac{T}{2}\sqrt{-\det(h)}h^{\eta\lambda}(\sigma)\frac{\partial g_{\mu\nu}}{\partial x^{\alpha}}\frac{\partial x^{\mu}}{\partial \sigma^{\eta}}\frac{\partial x^{\nu}}{\partial \sigma^{\lambda}} = 0$$
$$\frac{\partial}{\partial\sigma^{\rho}}\sqrt{-\det(h)} = \frac{1}{2}\sqrt{-\det(h)}h^{\gamma\beta}(\sigma)\frac{\partial h_{\gamma\beta}}{\partial \sigma^{\rho}} = \frac{1}{2}\sqrt{-\det(h)}\Gamma^{\gamma}_{\rho\gamma}$$
$$\begin{align}\biggr(-\frac{T}{4}\sqrt{-\det(h)}\Gamma^{\gamma}_{\rho \gamma}h^{\eta\rho}(\sigma)g_{\mu\alpha}(x(\sigma))\frac{\partial x^{\mu}}{\partial\sigma^{\eta}} + \frac{T}{2}\sqrt{-\det(h)}\frac{\partial h^{\eta\rho}}{\partial \sigma^{\rho}}g_{\mu\alpha}(x(\sigma))\frac{\partial x^{\mu}}{\partial \sigma^{\eta}} &+ \frac{T}{2}\sqrt{-\det(h)}h^{\eta\rho}(\sigma)g_{\mu\alpha}(x(\sigma))\frac{\partial}{\partial \sigma^{\rho}}\frac{\partial x^{\mu}}{\partial \sigma^{\eta}}\biggr) - \frac{T}{2}\sqrt{-\det(h)}h^{\eta\lambda}(\sigma)\frac{\partial g_{\mu\nu}}{\partial x^{\alpha}}\frac{\partial x^{\mu}}{\partial \sigma^{\eta}}\frac{\partial x^{\nu}}{\partial \sigma^{\lambda}} &= 0\end{align}$$
Does this make any sense?
Note: I am looking for a derivation without using Poincaré and Weyl invariance (see section "Symmetries of the Polyakov Action" and the comments). Details of the derivation can be found here, here and here, while this is for the Nambu-Goto action.
 A: $\newcommand{\d}{\mathrm{d}}\newcommand{\pd}{\partial}$Let's use worldsheet differential forms to make our life easier. So, for example,
$$\d X^\mu(\sigma) = \pd_\alpha X^\mu(\sigma)\; \d\sigma^\alpha.$$
I will also denote the fields by capital $X$; $\left\{\mu,\nu,\ldots\right\}$ will be spacetime indices and $\{\alpha,\beta,\ldots\}$ will be worldsheet indices (whenever necessary), so expressions like $\partial_\alpha$ will mean $\frac{\partial}{\partial\sigma^\alpha}$.
With these, the Polyakov action for the fields
$$X^\mu\in\operatorname{Maps}((\Sigma,h)\to(M,g))\cong \Omega^0(\Sigma;M)$$
can be rewritten as
$$S[X] = \frac{T}{2}\int_\Sigma g_{\mu\nu}(X)\ \d X^\mu\wedge \star\, \d X^\nu,$$
where $\star$ is the Hodge star operator on $\Sigma$. Note that from the point of view of the worldsheet, $\Sigma$, $g_{\mu\nu}$ is a scalar. It is then easy to just vary the action and get
\begin{align}
0\overset{!}{=} \delta S &= \frac{T}{2}\int_\Sigma \Big(\delta g_{\mu\nu}(X)\ \d X^\mu\wedge \star\, \d X^\nu + 2\,g_{\mu\nu}(X)\ \d\ \delta X^\mu\wedge \star\, \d X^\nu \Big) = \\ 
&= \frac{T}{2}\int_\Sigma \Big(2\, \Gamma^{(g)}_{\mu\nu\lambda}\ \d X^\mu\wedge \star\, \d X^\nu - 2\, \d\big(g_{\mu\lambda}(X)\ \star\, \d X^\mu\big) \Big)\ \delta X^\lambda,
\end{align}
where
$$\Gamma^{(g)}_{\mu\nu\lambda}:=\frac12\left(\frac{\partial g_{\mu\nu}(X)}{\partial X^\lambda}+\frac{\partial g_{\mu\lambda}(X)}{\partial X^\nu}-\frac{\partial g_{\lambda\nu}(X)}{\partial X^\mu}\right),$$
are the Christoffel symbols for $g$.
That's the EOM then:
$$ \d\big(g_{\mu\lambda}(X)\ \star\, \d X^\mu\big) -\Gamma^{(g)}_{\mu\nu\lambda}\ \d X^\mu\wedge \star\, \d X^\nu = 0.$$
You can reinstate the definition of the Hodge star and of $\d$ to get a more explicit presentation:
$$ \partial_\alpha\Big(\sqrt{-\det(h)}\, h^{\alpha\beta}\ g_{\mu\lambda}(X(\sigma))\;\partial_\beta X^\mu\Big) - \Gamma_{\mu\nu\lambda}^{(g)} \sqrt{-\det(h)}\, h^{\alpha\beta} \partial_{\alpha}X^\mu \partial_\beta X^\nu = 0.$$
Note that when $(M,g)=(\mathbb{R}^{D},\eta)$ this becomes the usual Polyakov EOM
$$ \partial_\alpha\Big(\sqrt{-\det(h)}\, h^{\alpha\beta}\;\partial_\beta X^\mu\Big) = 0.$$
The last one is simply ($\sqrt{-\det(h)}$ times) the Laplace-Beltrami operator on $(\Sigma,h)$ acting on $X^\lambda$. Using the formula for the variation of the determinant, one can show that this is equivalent to
$$ \sqrt{-\det(h)} \Big(h^{\alpha\beta} \partial_\alpha \partial_\beta - h^{\alpha\beta}{\Gamma^{(h)}}^\delta_{\alpha\beta}\partial_\delta\Big)X^\mu = 0,$$
where now ${\Gamma^{(h)}}^\cdot_{\cdot\cdot}$ are the worldsheet Christoffel symbols.
