I just tried to modify a code snippet for Neumann series to the Fredholm integral equation
$ f(x)= \sqrt{x} - \int^1_0 dy\; f(y) \sqrt{x y} $
which I read in another post: How to solve an integral equation by iteration method? with a slightly different integral, namely :
f[0]:=Sqrt[x];
f[n_]:=Sqrt[x]-Integrate[(f[n-1]/.x->y)*Sqrt[x*y],{y,0,1}];
data=Table[f[i],{i,0,10}];
The code works well but for my later application I need to solve those types of integrals numerically because the Kernel is too complicated to be solved by Integrate[]. Therefore I want to replace the integration by a numerical integration in the above code and naively tried
h[x_,0]:=Sqrt[x];
h[x_,n_]:=InterpolatingPolynomial[
Table[{x,Sqrt[x]-NIntegrate[h[y,n-1]*Sqrt[x*y],{y,0,1},AccuracyGoal->4]},{x,0,1,.1}], x];
data2=Table[h[x,i],{i,0,3}];
However even for a short series up to 3 iterations the computing time "explodes" and already exceeds the first code. How can I improve the code / computation?
I tried to replace the InterpolatingPolynomial by Interpolation[]
h[x_, 0] := Sqrt[x];
hint := Interpolation[Table[{x,Sqrt[x]-NIntegrate[h[y,n-1]*Sqrt[x*y],{y,0,1},AccuracyGoal->4]},{x,0,1,.1}]];
h[x_, n_] := hint[x];
data = Table[h[x, i], {i, 0, 3}];
but that code does not even compute the first iteration step h[x,1]