Let's start from the beginning:
The state of a quantum particle is represented by a vector $|\psi\rangle$ in the Hilbert space, observables are represented by Hermitian operators which eigenvalues represent the possible outcomes of a measurement and which eigenvectors represent the possible states after the measurement. One fundamental postulate is that the eigenvectors of any observable form a complete basis for the possible states of the particle. At my current level of understanding I think this is simply a postulate and cannot be demonstrated. Suppose we are dealing with an observable $R$ with only two possible results for a measurement: $+r,-r$, then we can write that the generic state of a quantum particle is: $$|\psi\rangle = a_+|+\rangle +a_-|-\rangle \ \ \ \ \ \ (1)$$ where $|+\rangle,|-\rangle$ are the eigenvectors of $R$. The probabilities of one outcome or the other are respectively: $|a_+|^1,|a_-|^2$, this is also a postulate.
But what to do if our observable is, for example, the position, with an infinite number of possible measurement outcomes? We start by noticing that we can rewrite (1) as: $$|\psi\rangle = \langle +|\psi\rangle|+\rangle +\langle -|\psi\rangle|-\rangle \ \ \ \ \ \ (2)$$ this is just a mathematical property of the Hilbert space. We can then think to have an infinite number of basis vector with an infinite number of coefficients: $$\psi=\langle x_1 | \psi\rangle|x_1\rangle+\langle x_2|\psi\rangle|x_2\rangle+.....$$ this reasoning is not really rigorous, but in the limit of dense and infinite eigenstates allows us to think about the collection of coefficients as a function: $$\psi(x)=\langle x |\psi\rangle$$ then to find the probability of measuring $x$ in an interval $(a,b)$ we can think to sum all the probabilities of $x$ in that region, but since we are in the mentioned limit we integrate instead of summing: $$P(a<x<b)=\int_a^b|\psi(x)|^2dx$$ From what I currently understand this is what a wave function is!
Long story short: for me the wavefunction of a particle is a complex valued function that when integrated in its variable/s gives the probability of finding the particle under the interval of integration.
But I started to doubt my understanding, I think that the wave function may be more versatile than this. To explain why lets take an exercise as an example:
Given the wave function: $$\psi(r,\theta,\phi)=Ae^{-br}(1+2br\sin{\theta}\sin{\phi})$$ Find the possible outcomes of a measurement of $L^2,L_z$ with relative probabilities.
I strongly suspect that what I am about to say will appear ridiculous to a seasoned expert on this topic, but unfortunately I am not an expert, so hear me out:
Given my understanding of what a wave function is this question does not make sense! If you give me a wave function in the variables $r,\theta,\phi$ I can then give to you the probabilities for the particle to be in an area of $3D$ space. Because this is what a wave function allows us to do by definition, this is the information that it's carrying by construction! But how and why can I infer information about $L^2$ and $L_z$ with the wave function only?
This exercise of course has a standard procedure of resolution: write the wavefunction in terms of the eigenfunctions of $L^2,L_z$, so the Spherical Armonics $Y_{l,m}(\theta,\phi)$; if we do this we get: $$\psi=Ae^{-br}\sqrt{4\pi}\left[Y_{0,0}+\sqrt{\frac{2}{3}}ibr(Y_{1,1}+Y_{1,-1})\right]$$ and then we can somehow find the probabilities of $l=0,m=0$, ecc. by integration. Problem is I don't understand how or why this method works! I mean: the rewriting in terms of $Y(\theta,\phi)$ is mathematically ok, but why is it useful? The wave function, written in one way or the other, should still carry only the information about the probability of finding the particle in a region of space. How can we extract the information we need and why?
Is my definition/understanding of the wavefunction wrong? Is it more versatile than I thought? How?