I have a recurrence relation as follows:
$$ d(n+2) = -(n+2)^2d(n) - (2n+5)d(n+1)$$
Setting n=0 and generating a few coefficients gives
$ d(0) = a$
$ d(1) = b $
$ d(2) = - 4a - 5b $
$ d(3) = 28a + 26b $
$ d(4) = - 188a - 154b $
$ d(5) = 1368a + 1044b $
$ d(6) = - 11016a - 8028b $
It turns out the coefficients of a and b are both listed in the OEIS (a, b). Both mention being a Harmonic number times factorial, but the connection doesn't seem very rigorous (e.g., "If we eliminate the first term and set a(0)=3, then a(0)..a(12) = (-1)^(n+1)n!(2*H(n)-3)..where H(n) is the n'th harmonic number").
Given the initial conditions d(0)=a, d(1)=b Wolfram alpha gives
$$ d(n) = -(-1)^n \Gamma(n+2) ((2 a+b) H_{n+1}-3 a-b) $$
I'm trying to get from the recurrence equation to the slightly more accessible form with Harmonic numbers, but I don't know how.
This was answered in Marko Petkovsek's Thesis (who also implemented the algorithm in mathematica). Presumably (but not certainly), a more expository account is given in Petkovsek/Zeilberger's A=B (which is available for free).
EDIT The Petkovsek paper discusses higher order linear difference equations too.