rsolve.sage

def rsolve(coeff, values, rhs=0, nonhom_sol=0):
    R.<x> = CC[]
    f = 0
    for i, c in enumerate(coeff):
        f += c * x**i
    h = 0
    n = var('n')

    cs = []
    cnt = 0
    for r, m in f.roots(SR):
        for i in range(m):
            cnt += 1
            c = var('c%d' % cnt)
            cs += [c]
            h += c * n**i * r**n

    g = h + nonhom_sol
    #print g
    eqs = []
    for arg, val in values:
        # there is a bug in solve: solve([x==1], [x], solution_dict=True) fails
        # so just add equation twice
        eqs.append(g(n=arg) == val)
        eqs.append(g(n=arg) == val)
    sol = solve(eqs, cs, solution_dict=True)[0]
    for c in cs:
        g = g(**{str(c):sol[c]})

    return g

n=var('n')
g = rsolve(
    [-1,0,1,1,1,-1,-1,-1,0,1],
    [(0,0), (1,0), (2,0), (3,0), (4,1), (5,1), (6,2), (7,3), (8,5)],
    1,
    1/144*n^3
    )
print g
for i in range(100):
    print i, g(n=i).expand()