So the way I understand the question is that it asks us at what r (distance from nucleus) and at what angles is the probability going to be maximum. The way I see probability of finding at a point in a meaningful sense is, just like in 1D probability of finding particle at a point is taken as probability it is in dx region beside it, we find probability of finding it in a small volume near the point of interest. Now the probability density of 1s orbital times dV in polar coordinates will end up containing a sine term (since small volume in polar cooditnates contains sine of angle with z axis). Thus the probability is dependent on not just r but also angle with z axis, which isn't expected for such a spherically symmetric orbital. I am confused here.
Of course one could say that since this is spherically symmetric, probability depends only on r and thus use only radial probability distribution and maximise it. But why is the above way which seems more general giving an apparently different answer.
Sorry if the above part isn't clear or is messy. Just to clarify things, I am not asking a "how to find" question. One way to look at my question is that if instead of maximising radial probability distribution, if I evaluate probabilities in small volume around points and maximise it, my answer should match with the one I get from maximising radial probability distribution. Consequently, I should get all points on sphere with Bohr radius but on maximising probability in small volume around a point, I get points at Bohr radius and on XY plane only (since for maxima, sine must be 1 meaning perpendicular to z axis).