[EDIT]
Hold on, I was confused by theta/phi ambiguity and your component ordering.
You have the correct spherical vector but in an unfamiliar order.
Since I typed all of this I will not leave it in case anybody else is helped by it.
What is p=rsinθ? the value of function f or the radial unit variable in spherical because you used r not p in the formula. I use ρ,φ,θ in spherical to distinguish from r, θ, z
Almost right. pretty close. Cylindrical looks good. But Spherical is off in the two angular variables. You might beak it up into two operations. Laplacian is basically the divergence of the gradient of a scalar field. In cylindrical the azimuthal differential , rdθ , depends on r and θ because to move in the θ direction length is an arc depending on radius. In spherical your surface differential depends on r, Φ, and θ. to get an area on the surface of a sphere, coming down from north pole your distance is an arc with differential length rsinθ just like cylindrical turned vertically but then the other dimension is rsinΦ or r^2sinθsinΦ. That is the differential unit surface.
It looks like you have f a function of r and θ. If this is the case the you need to keep them in and do the derivatives instead on pulling them out of the del operator like they were constants. if your function is rsinθ you have 1/r^2 d/dr( r^2 d/dr( rsinθ)) then 1/r^2 d/dr( r^2 sinθ) last 2sinθ/r. this is the final r hat component. In the formula you posted the derivatives were shown as two nested derivatives instead of the double derivative d^2/dr^2 for this reason. The same for the zenith variable (down from z axis, some use theta some use phi). For the azimuth variable (around the equator) this is worked out.
