(* This is going to solve the discretised integral equation for Exercise 3.8 *)

(* First set up and solve the linear system *)

A[n_,m_,d_,e_]:=ArcSinh[d(m-n+1/2)/e/S[(n+1/2)d]]-ArcSinh[d(m-n-1/2)/e/S[(n+1/2)d]]

prob[e_,max_]:=Table[f[n]-Sum[A[n,m,1/max,e]f[m],{m,-max,max-1}]==4Pi,{n,-max,max-1}]

sol[e_,max_]:=Table[{(n+1/2)/max,f[n]},{n,-max,max-1}]/.Evaluate[Solve[prob[e,max],Table[f[n],{n,-max,max-1}]][[1]]]

(* Print solution and compare with leading-order approximation
f(z) ∼ -((4Pi)/(2log(2/ϵ)-1)) *)

printsol[e_,max_]:=Show[Plot[-4Pi/(2Log[2/e]-1),{z,-1,1},PlotStyle->Darker[Red],PlotRange->Full],ListPlot[Evaluate[sol[e,max]]],Plot[-4Pi/(2Log[2/e]-1),{z,-1,1},PlotStyle->Darker[Red]],AxesLabel->{z,f[z]}]

(* First example: ellipsoid *)

S[z_]:=Sqrt[1-z^2]

printsol[0.02,50]

(* Pretty good. Approximation should be exact in this case *)

(* Another example *)

S[z_]:=Cos[Pi z/2]

printsol[0.02,50]

(* Approximation not so close but still gives right order of magnitude *)