Consider a Schwarzchild black hole of mass M=15Ms where Ms is the mass of the sun

I'll explain how to derive the "free fall formula" for radial coordinate as a function of proper time for a free falling object in the radial direction.

The equation of motion follows from the Lagrangian.

L = 1/2 g_{mu,nu}x^{mu}-dot x^{nu}-dot

where the dot means the derivative w.r.t. proper time.

The action is Integral over tau of L dtau and the Euler Lagrange equations are thus:

d/dtau [dL/dx^{mu}-dot] - dL/dx^{mu} = 0

You can read-off the g_{mu, nu} from the line element:

ds^2 = c^2(1-2m/r)dt^2 - (1-2m/r)^(-1)dr^2 - r^2 d theta^2 - r^2sin^2(theta)d phi^2

where m = MG/c^2. So, L is given by:

L = c^2(1-2m/r)t-dot^2 - (1-2m/r)^(-1)r-dot^2 - r^2 theta-dot^2 - r^2sin^2(theta) phi-dot^2

Because ds^2 = c^2 dtau^2 you also have the equation:

c^2(1-2m/r)t-dot^2 - (1-2m/r)^(-1)r-dot^2 - r^2 d theta-dot^2 - r^2sin^2(theta)d phi-dot^2 =c^2

If we are free falling in the radial direction, then theta-dot = phi-dot = 0, so this equation simplifies to:

c^2(1-2m/r)t-dot^2 - (1-2m/r)^(-1)r-dot^2 = c^2 (1)

We now only need to eliminate t-dot from this equation. So we write down the Euler-Lagrange equation for mu = 0 (i.e. the time coordinate):

d/dtau dL/dt-dot - dL/dt = 0

L does not depend on t, and we find that:

d/dtau[(1-2m/r)t-dot] =0 -->

(1-2m/r) t-dot = k (2)

where k is some constant.

Using (2) you can eliminate t-dot in (1):

1/2 r-dot^2 - mc^2/r = 1/2 c^2 (k^(2) - 1) (3)

If the particle is released at coordinate r0, then r-dot must be zero for r = r0 and we can eliminate k in the above equation:

1/2 r-dot^2 - mc^2/r = - mc^2/r0

m = MG/c^2, so we can rewrite this as:

1/2 r-dot^2 = MG[1/r - 1/r0] (4)

From this equation you find the proper time needed to fall from r0 to r1:

tau_{r0,r1} = 1/sqrt[2MG]Integral from r1 to r0 of dr sqrt[r0*r/(r0-r)]

You can easily calculate this integral by substituting r = r0sin^(2)(u).

If you substitute the radial coordinate back in the result you should obtain:

tau{r0,r1} = r0^(3/2)/(sqrt[2MG]) {pi/2 - arcsin[sqrt(r1/r0)] + sqrt[(r1/r0) (1-r1/r0)]}