How to make this algorithm stable

Question

I have to compute the iterative version of this formula:

f(i)=integral ( x^i/(4x+1) ) from 0 to 1

Using these formulas:

f(0)=ln(5)/4;
f(i)=1/(i+1) - 4*f(i+1);

I have tried tried the following: I calculate integral ( x^100/(4x+1) ) from 0 to 1, and store the result.Then I calculate f(i) starting from this result, using the iterative version.
But I get wrong result because the error is too big.
The error is acceptable only for i<=25.
I would like to know, why this algorithm is not stable and if there's a solution to compute the result starting from i=100 or higher.

This is the code:

function y=Integral(i,max)

if i==0
    y=1/4*log(5);
elseif i==max
    y=0.0;
else
    y=1/(i+1)-4*Integral(i+1,max);
end


end

With this function I never get an exact value because the error accumulated is too high.I get a close value (but even 3 or 4 times higher, so not acceptable) if I use i=15 and max=18.I need the stable version of this formula.

Answer 1

This recursive function should do the job without the need to store partial results on the way to 100:

function y = Integral(i)

if i==0
    y=log(5)/4;
else
    y = (-Integral(i-1) + 1/i)/4;
end
end

For recursion to work you need to start from i=100 and then call the function with i-1 until it reaches i=0.

Integral(100) will give the final answer without the need to store partial results.