(* First define a useful function *)
Useful[k_]:= k^2 (k + 1)^2 (2 k^2 + 2 k -  1);
(* Now define the denominator of P[n] *)
Q = Useful[M]

(* Now define P[n] itself *)
P[n_]:= Useful[n]/Q;
(* Write out P[n] *)
P[n]

(* Write out P[n-1] *)
P[n-1]

(* Now find P[n] - P[n-1] *)
P[n] - P[n-1]

(* What a mess! Simplify it *)
Simplify[P[n] - P[n-1]]

(* Simplify a bit more by factoring *)
Factor[Simplify[P[n] - P[n-1]]]

(* Enter (39) *)
approximation = t^6 ((Useful[n-1]/12)/n^6 + 1/(2n))

Simplify[%]

(* As n -> Infinity, that thing -> t^6/6 *)