Hey, I would like know how to find the field amplitude with equation 26 shown in figure 6 using the listings below:

 d = 1; lam = d; a = 25*d; b = 100*d ; x2 = 0.5 a; k = 2 \[Pi]/d
    u[0][x_] := 1; v[0][x_] := 0
    Do[{u[q] = 
        Interpolation[
         Re[c*NIntegrate[(u[q - 1][x1] + 
               I*v[q - 1][x1])*(Exp[-I*k*a*Sqrt[b^2 + (x1 - x2)^2]]/
               Sqrt[Sqrt[b^2 + (x1 - x2)^2]])*(1 + 
               b/Sqrt[b^2 + (x1 - x2)^2]), {x1, -1, 1}]]], 
       v[q] = Interpolation[
         Im[c*NIntegrate[(u[q - 1][x1] + 
               I*v[q - 1][x1])*(Exp[-I*k*a*Sqrt[b^2 + (x1 - x2)^2]]/
               Sqrt[Sqrt[b^2 + (x1 - x2)^2]])*(1 + 
               b/Sqrt[b^2 + (x1 - x2)^2]), {x1, -1, 1}]]]}, {q, 1, 10}];
    amp = Interpolation[
       Table[{q, (Abs[u[q][.5] + I*v[q][.5]])^2}, {q, 0, 10}]];
    {Show[{ListPlot[
        Table[{q, (Abs[u[q][.5] + I*v[q][.5]])^2}, {q, 0, 10}], 
        PlotRange -> All, Frame -> True, 
        FrameLabel -> {"Number of Transits", "Amplitude"}], 
       Plot[amp[q], {q, 0, 10}, PlotRange -> All]}], 
     Plot[Table[(Abs[u[q][x] + I*v[q][x]])^2, {q, 1, 10, 1}], {x, -1, 1}, 
      ColorFunction -> Hue, Frame -> True, 
      FrameLabel -> {"x", "Amplitude"}, Axes -> False]}

Cordially.

Attachments: