The probability distribution function is the square of the absolute value of the wavefunction. As you approach the classical limit the wavefunction will oscillate faster and faster. To get a well defined classical probability distribution, you define a coarse-grained average by integrating over a small length interval epsilon before you take the classical limit. This is not just a mathematical trick. If you want to measure the position of the particle you must use some device that has a finite length. So, you will always have to deal with the average of the probability distribution over some finite length interval.If you take the classical limit then this coarse grained average will be smooth function because you've averaged out the infinite fast oscillations. You can then take the limit epsilon --> 0 to zero without problems. That function will still be smooth and this is the classical probability distribution P(x).Let's start with deriving the energy eigenfunctions.
Outside the square well the wavefunction must be zero. Because the wavefunction is continuous the wavefunction is zero at the boundaries. Inside the square well the solutions are linear combinations of exp(i k x) and exp(-i k x) where the energy is given byE = hbar^{2} k^{2}/(2m)If one of the boundaries is at x=0 then, working in terms of trigonometric functions, you must discard the cosine term because that term is 1 at x = 0. So, you are left with a sinus term:A sin(kx) If the other boundary is at x=L, then demanding that the wavefunction is zero there yields:
k L = n pi The wavefunction is thus: A sin(npi x/L) The square of this function integrated from zero to L must be 1. This implies that: A = sqrt[2/L] the normalized wavefunctions are thus given by: psi_{n,L}(x)=sqrt[2/L]sin[n pi x/L] (1)

and the energy is: E_{n,a}= hbar^{2}pi^{2}n^{2}/[2 m L^{2}] (2) For a classical objects, the mass and the energy will take values that are not very small when expressed in SI units. This means that the quantum number n must be very large. This means that the wavefunction (1) must oscillate very rapidly. You can also take the classical limit by sending h-bar to zero while keeping the energy constant. This has the same effect: n goes to infinity.

The probability distribution is the square of psi_{n,L}(x):

P_{n,L}(x) = 2/L sin^{2}[n pi x/L] (3)

The coarse-grained average is the integral of this function from x to x + epsilon divided by epsilon. You don't need to work this out explicitly for general values of n as we only need it for the limit n to infinity. In this limit the average of the sin^{2} term becomes 1/2. To see this note that if you choose: n>> L/epsilon (4) then there will be a large number of oscillations inside a length element of epsilon. The average of sin^{2} over an integer number of periods is exactly 1/2. The maximum deviation from 1/2 is thus less than the maximum value the sin^{2} can attain times the period of the function divided by epsilon. This is of order L/(n epsilon) and from (4) you see that this must be much smaller than 1. So, we can safely replace the Sin^{2} term in (3) by 1/2 and obtain the probability distribution:

P(x) = 1/L (5)

This no longer depends on epsilon and the limit epsilon --> 0 is thus trivial.

The average of x is the integral of xp(x) from zero to L and this is L/2.

The average of x^2 is the integral of x^{2}p(x) from zero to L and this is L^{2}/3.

If you need the averages w.r.t. the center of the square well, then you can use that the average of x is zero w.r.t. the mid point (the average of x calculated above is L/2 which is the mid point, but you can also reason as follows. Because the probability distribution is symmetric w.r.t. the mid point and the position relative to the mid point is an odd function, the average must be zero). Now, in general you have:

<x-<x>>^{2} = <x^{2}> - <x>^{2} (6)

Since <x> is the mid point, x - <x> is the position relative to the mid point. Eq. (6) says that the average of the square distance measured from the mid point is

<x^{2}> - <x>^{2}

inserting <x^{2}> = L^{2}/3 and

<x> = L/2 in here gives:

<x-<x>>^{2} = L^{2}/12

So, if we measure the position relative to the center, the average position is zero and the average of the square of the position is L^{2}/12