Working with the central binomial coefficients squared.
Requirements: ore_algebra package
# Define the ring of shift operators to encode P-recurrence
from ore_algebra import *
Ind.<n> = PolynomialRing(QQ); Shift.<Sn> = OreAlgebra(Ind)
# Define the P-recurrence satisfied by the sequence
rec = (n+1)^2*Sn-4*(2*n+1)^2
rec
(n^2 + 2*n + 1)*Sn - 16*n^2 - 16*n - 4
# Compute list of sequence terms using recurrence (using that initial term is 1)
rec.to_list([1],10)
[1, 4, 36, 400, 4900, 63504, 853776, 11778624, 165636900, 2363904400]
# Verify we encode the central binomial coefficients squared
[binomial(2*n,n)^2 for n in range(10)]
[1, 4, 36, 400, 4900, 63504, 853776, 11778624, 165636900, 2363904400]
# Define ring of linear differential operators
Pols.<z> = PolynomialRing(QQ); Diff.<Dz> = OreAlgebra(Pols)
# Compute a D-finite equation satisfied by the GF of the central binomial coefficients
diff = rec.to_D(Diff)
diff
(-16*z^3 + z^2)*Dz^3 + (-80*z^2 + 3*z)*Dz^2 + (-68*z + 1)*Dz - 4
# Note the textbook example has a differential equation of smaller order
textbook_diff = z*(1-16*z)*Dz^2 + (1-32*z)*Dz - 4
textbook_diff
(-16*z^2 + z)*Dz^2 + (-32*z + 1)*Dz - 4
# The differential operator calculated here is a left-multiple of the textbook operator
diff.quo_rem(textbook_diff)
(z*Dz + 1, 0)
# In other words, diff = (z*Dz + 1)*textbook_diff.
# Because diff has a solution which begins with the same initial terms
# as the generating function, the GF is also annihilated by textbook_diff
diff - (z*Dz + 1)*textbook_diff
0