The problem says
Show that when Laplace's equation $$ \frac{\partial ^2 u}{\partial x^2}+\frac{\partial ^2 u}{\partial y^2}+\frac{\partial ^2 u}{\partial z^2}=0 $$ is written in cylindrical coordinates, it becomes $$ \frac{\partial ^2 u}{\partial r^2}+ \frac{1}{r} \frac{\partial u}{\partial r} + \frac{1}{r^2} \frac{\partial ^2 u}{\partial \theta^2} + \frac{\partial^2 u}{\partial z^2}=0 $$
I've been stuck on this problem for hours and I can't find a relatively elegant way of solving this. I started by restating the problem into showing that $$ \frac{\partial ^2 u}{\partial x^2}+\frac{\partial ^2 u}{\partial y^2}+\frac{\partial ^2 u}{\partial z^2} = \frac{\partial ^2 u}{\partial r^2}+ \frac{1}{r} \frac{\partial u}{\partial r} + \frac{1}{r^2} \frac{\partial ^2 u}{\partial \theta^2} + \frac{\partial^2 u}{\partial z^2} $$ And the z coordinate just cancels out, but I hate these types of problems because they evolve into an unreadable and confusing mess with a couple of applications of the chain rule. Also, this problem is placed after a unit on multiple integrals, but I can't connect how multiple integrals are related to this problem. The closest thing that I can connect to is cylindrical coordinates (who would have thought) and the Jacobian (but I don't think that it will be very useful here). Can anybody give me a hint (please don't just solve the entire problem and leave the last step for me as a "hint", I just want a nudge in the right direction, so a vague question would be great) for a practical way of solving this? Or is "brute force" (applying the chain rule and all that) the only way through this?
Go back to basic definitions. Over a Riemannian manifold, the Laplacian is the dual with respect to the bilinear integral form on scalar products of gradients for two functions $f,g$, being elements of a kernel of smooth functions with compact support away from boundaries.
From such a core the algebraic form is deduced, that is taken to be the definition of the Laplacian for more general scalar functional object classes.
$\int dV \ ( f \ \Delta g ) = - \int dV\ ( \nabla \ f \cdot \nabla \ g ) $
Since for functions, partial and covariant deriviatives coincide, the algebraic scheme works quite generally for all forms of n-forms.
In cylindrical coordinates the scale factors are $h_r=1,\ h_z=1,\ h_\phi= r$, the volume is the product in orthogonal coordinates $dV = r\ dr \ dz\ d\phi $
$\int \ r \ dr\ dz\ d\phi\ ( f \ \Delta g ) = - \int r\ dr \ dz \ d\phi \ (\partial_r f, \ \partial_z f, \ r^{-1} \partial_\phi f ) \ \cdot \ (\partial_r g, \ \partial_z g, \ r^{-1} \partial_\phi g ) = \int r \ dr \ dz \ d\phi \ \bigl( f\ \ (r^{-1} \ ((\partial_r,\ \partial_z, \ r^{-1} \partial_\phi))) \ \cdot \ r \ ((\partial_r , \ \partial_z,\ r^{-1} \partial_\phi))) g \bigr) $
$\Delta g = r^{-1}\partial_r(\ r \ \partial_r\ g) + \partial_{zz}\ g + r^{-2} \partial_{\phi \phi}\ g $
Obviously the general formula in an orthogonal system of coordinates is
$\Delta = \sum_i \frac{1}{\prod_k h_k} \partial_i (\frac{\prod_m h_m}{h_i^2} \partial_i) $
The interplay of scale factors and volume elements in dualization generally reflects the volume expansion in the different directions, that count fundamentally for the question of the divergence as the source density per standard euclidean unit volume.
Becuase the algebra of derivatives of the coefficients in curvilinear coordinates in dimensions $>2$ is rather chaotic, its always better to use the tables.