My problem is that I want to solve the integral equation $$ f\left(t\right) = K\left\lbrace e^{\int_{t}^{T}\left(\alpha + \beta f(s) + \sigma f^2(s)\right)ds} + \gamma e^{-\frac{1}{2}\int_{t}^{T} \sigma f^2(s) ds} - \gamma \right\rbrace, $$ where $K,\alpha,\beta,\sigma,\gamma >0$ are all constants. One can show that if we construct a sequence $\left\lbrace f_n(t) \right\rbrace_{n\geq 0}$ with $f_0=1$ and $$ f_{n+1}(t) = K\left\lbrace e^{\int_{t}^{T}\left(\alpha + \beta f_n(s) + \sigma f_n^2(s)\right)ds} + \gamma e^{-\frac{1}{2}\int_{t}^{T} \sigma f_n^2(s) ds} - \gamma \right\rbrace, $$ it converges, $f_n \to f$ as $n\to \infty$. That is, we can solve the initially stated equation numerically by iteration starting from $f_0=1$.
My question is, how would you perform these iterations? I tried writing a Matlab algorithm, and it can't handle more than 3 iterations. Mathematica stops at 4.