How do I use RSolve to solve a system of recurrence relations?

I am trying to solve a system of recurrence relations as follows.

RSolve[{a[p] == 1 + (2/n) b[p - 1] + ((p - 1)/n) a[p - 1],
b[p] == 1 + (p/n) b[p - 1],
a[1] == 1 + 2/n,
b[0] == 1},
{a[p], b[p]}, p]


Having done this, I would like to set $p=n-1$ and then plot the resulting functions $a[n-1]$ and $b[n-1]$ for $n=1, ..., 1000$.

I am failing at the first hurdle as Mathematica version 8 just returns the RSolve line when I press shift-enter.

Update. Following a request for an example using Maple 16. I do the following.

rsolve({a(1) = 1+2/n, a(p) = 1+(p-1)*a(p-1)/n+2*b(p-1)/n, b(0) = 1, b(p) = 1+p*b(p-1)/n}, {a(p), b(p)});
subs(p = n-1, %);


I then copy the result for a(n-1) and plot it

plot([seq([n, (1/n)^(n-1)*n*GAMMA(n-1)*(Sum((n+2*(1/n)^p1*exp(n)*GAMMA(p1+1, n))/(n^2*(1/n)^(p1+1)*GAMMA(p1+1)), p1 = 1 .. n-2)+1+2/n)], n = 2 .. 200)], style = point);


This gives

-
First issue is the syntax : RSolve with capital S. –  b.gatessucks Jan 1 '13 at 19:30
Strongly related mathematica.stackexchange.com/q/16919/193 –  belisarius Jan 1 '13 at 19:41
Sorry... copy and paste error. I have the same problem with RSolve with a capital S. Is there a way to restart so that there are no existing function/variable settings or is there some other reason it doesn't appear to do anything? –  lip1 Jan 1 '13 at 19:48
solving for b[p] seperately produces a solution for b : solb = RSolve[{b[p] == 1 + (p/n)*b[p - 1], b[0] == 1}, b, p][[1, 1]]. Substituting this in the other equation seems to be more complicated though sola = RSolve[{a[p] == 1 + (2/n)*b[p - 1] + ((p - 1)/n)*a[p - 1], a[1] == 1 + 2/n} /. solb, a, p] –  Thies Heidecke Jan 1 '13 at 19:53

In two steps, since they can be decoupled:

h1 = RSolve[{
b[p] == 1 + (p/n) b[p - 1],
b[0] == 1}, {b[p]}, p]
h2 = RSolve[{
a[p] == Evaluate[b[p] /. h1[[1]]] + ((p - 1)/n) a[p - 1],
a[1] == 1 + 2/n}, {a[p]}, p]

DiscretePlot[{b[i] /. Evaluate[h1 /. p -> n - 1 /. n -> i + 1], Sqrt[i]}, {i, 300}]
DiscretePlot[ a[i] /. Evaluate[h2 /. p -> n - 1 /. n -> i + 1], {i, 300}]


Comparing b[i]with Sqrt[i]

And this is a[i]

Edit

Thanks. How can I now do the substitution and plotting? If there isn't a bug then it should look like $sqrt(n)$

here is a reasonable approximation to b[i] using Sqrt[]:

y[i_] := Sqrt@(1.78081659706051 + 1.57076044099885*i) - 0.332147990949542;
DiscretePlot[Evaluate[b[p] /. h1[[1]] /. p -> n - 1 /. n -> i + 1] - y[i], {i,10000},
PlotRange -> All]


Or, if you want an exact result:

b[n] -> E^n n ExpIntegralE[1 - n, n]

-
Thanks. How can I now do the substitution and plotting? If there isn't a bug then it should look like $\sqrt{n}$. –  lip1 Jan 1 '13 at 19:56
Thanks but isn't it possible to use the answer to RSolve[{ a[p] == Evaluate[b[p] /. h[[1]]] + ((p - 1)/n) a[p - 1], a[1] == 1 + 2/n}, {a[p]}, p] without having to retype the answer for a[n]? –  lip1 Jan 1 '13 at 20:04
@lip1 Look at it now –  belisarius Jan 1 '13 at 20:13
Thanks very much. –  lip1 Jan 1 '13 at 22:30
@belisarius Don't you think Mathematica should be able to solve it the way you did ? –  b.gatessucks Jan 2 '13 at 7:59