It's a standard example, but I thought somebody might enjoy the Maxima
way.
We need nfloat, so load hypergeometric.
(%i1) load(hypergeometric)$
(%i2) i(n) := if n = 0 then 1-1/%e else n * i(n-1) - 1/%e$
(%i17) [float(i(42)),
float(expand(i(42))),
block([fpprec : 52], bfloat(i(42))), block([fpprec : 56],
bfloat(i(42))),
nfloat(i(42),[],20)];
(%o17) [-4.7467928326479132*10^34,
0.0,
2.732067241502428018740764014347935258990086011897504b-2,
8.7536981302551482132924825234864659298743235134456890016b-3,
8.7541921419809776312b-3]
Solving the recursion for I(42) is identical to solving a linear
system with coefficient matrix
(%i19) m : genmatrix(lambda([i,j], if i-j = 0 then 1 else if i - j = 1
then -i else 0), 42,42)$
And the condition number of m is big--and I mean big
(%i22) mat_norm(m,inf) * mat_norm(m^^-1, inf);
(%o22) 103810448683368212798606246049077583803159875197431095
Barton