Waiting some time for a straightforward numerical answer, here my attempt, which assumes ` \[Beta] = 1` (without loss of generality) and predefined values `c0==0.2,T0==5` : c0 = 2/10; T = 5; gip[eps_] :=Module[{x, t}, Interpolation[Table[{t,eps[t]^2 + c0 + Exp[-t] NIntegrate[eps[s] Exp[s], {s, 0, t}] -Exp[-t] (1 + eps[t]) NIntegrate[eps[s]^2 Exp[s], {s, 0, t}]}, {t, Subdivide[0, T, 25]}]]] Function `gip` get's a pure function as input argument and returns an InterpolationFunction, which might be used iteratively. With appropriate starting value `eps[0]== 1/2 (1 - Sqrt[1 - 4 c0])` (second solution branch `eps[0]== 1/2 (1 + Sqrt[1 - 4 c0]) omitted` ) it follows solm = NestList[gip, 1/2 (1 - Sqrt[1 - 4 c0]) &, 15] ; Plot[Evaluate[Through[solm [t]]], {t, 0, T}, PlotRange -> {0, All}] [![enter image description here][1]][1] The solution doesn't match completely @cesaro's one (shows superimposed oscillations). Perhaps it shows a first step for a completely numerical solution. To exclude the influence of simple `Interpolation` I also tried a numerical solution using `NDSolveValue` (substitution `i1[t]==Integrate[eps[s]Exp[s],{s,0,t}]` and `i2[t]==Integrate[eps[s]^2 Exp[s],{s,0,t}]`) : ff = Function[{t}, NDSolveValue[{i1'[t] == #[t] Exp[t], i1[0] == 0,i2'[t] == #[t]^2 Exp[t], i2[0] == 0},#[t]^2 + c0 + Exp[-t] (i1[t] - (1 + #[t]) i2[t]) , {t, 0,T}, DependentVariables -> {i1, i2},Method -> {Automatic ,"DiscontinuityProcessing" -> False } ,AccuracyGoal -> 10 ]] & ; which gives the same result(I used only 7 iterations because of increased evaluation time): sol = NestList[ff, 1/2 (1 - Sqrt[1 - 4 c0]) &, 7]; // AbsoluteTiming Plot[Evaluate[Through[sol[t]]], {t, 0, T}, PlotRange -> {0, All}] [![enter image description here][2]][2] **conclusion** The two numerical solutions, obtained with two independent methods, match very well. I'm quite convinced that the solutions describe the iteration of the given integral equation. Hints for improvement are welcome! [1]: https://i.sstatic.net/cRK1H.png [2]: https://i.sstatic.net/8qZRh.png