\(\(\( (*\ 
      This\ Mathematica\ notebook\ computes\ the\ asymptotic\ number\ of\ \ pi' 
        shapes\ which\ are\ compatible\ with\ an\ RNA\ sequence\ of\ length\ n, \ under\ the\ assumption\ that\ any\ base\ can\ basepair\ with\ any\ other\ base\ and\ that\ there\ are\ at\ least\ 3\ unpaired\ bases\ in\ a\ hairpin\ \(\(loop\)\(.\)\)\ *) \)\(
\)\(
\)\(\(Clear["*"];\)

    eqn = {U == z + 1, \ 
        S == U*S*T*z^2 + U, 
T == 
          U^2*T^2*S*z^4 + 2*z^3*T + z^4*T + z^3, 
S0 == 
          z*S0 + S}

    solset = Simplify[Solve[eqn, S0, {U, S, T}]]



    \)\)\)

One of these 2 is the generating function.  Setting z=0, we can quickly see that first one blows up as z->0 so must be second one.

\(\(\(sol = S0 /. Extract[solset, 2]\)\(
\)
    \)\)

Singularity occurs either where denominator is zero (0 or ±1, but we know analytic at z=0 as it is a generating function.) or where square root is zero.  Let's find where square root is zero.

poly = 
      1 - 4 (z^3) - 2 (z^4) - 2 (z^5) + 2 (z^6) - 7 (z^8) - 6 (z^9) - (z^10) + 
        2 (z^11) + (z^12)

roots = NSolve[poly == 0,  z]

r = z /. Part[roots,  9,  1]

This is our (unique) dominant singularity.  (It is smaller in magnitude than the singularities z=±1 caused by the denominator.)

The next few commands are to simply divide poly by (1-z/r) to get poly2 and then check that this worked.

poly2a = poly/(1 - z/r)

Series[poly2a,  {z,  0,  14}]

poly2 = 
      1 + 1.8077595903919128` z + 3.2679947366539364` (z^2) + 
        1.907748826536447` (z^3) + 1.4487512372301783` (z^4) + 
        0.6189939431950044` (z^5) + 3.1189922372052834` (z^6) + 
        5.638388129165776` (z^7) + 3.1928502148513473` (z^8) - 
        0.22809440341762866` (z^9) - 1.4123398452929337` (z^10) - 
        0.5531709002208345` (z^11)

Expand[poly2*(1 - z/r)]

This is close to poly, so our polynomial division worked.

poly

left = (1/(2 (z^2) (-1 + (z^2))))*Sqrt[poly2]

This is the portion next to the singularity, (1-z/r)^(1/2).

K = left /. z -> r

K2 = K/Gamma[-1/2]

1/r

So V_n ~ 1.27613*(1.80776)^n*n^(-3/2)

polySol = Series[sol,  {z,  0,  2000}];

We can use the above series to check our asymptotics.  Let's pull out the term in front of the z^500.

coeffs = CoefficientList[polySol,  z];

act500 = Part[coeffs,  501]

est500 = K2*(1/r)^500*500^(-3/2)

N[act500/est500]

This shows our asymptotics to be close at n=500.  By the same method we can check at other n, or write a program to compare for all n.