in these functions deltas at +-a; in web delta at +-a/2; plotting uses web version 1/a or a on x axis
fs[a_]=-1/2-Exp[-3 a](Exp[-a]+Exp[+a])/(1+2(1/2+a)Exp[-2 a])
fa[a_]=-1/2-Exp[-3 a](Exp[-a]-Exp[+a])/(1-2(1/2+a)Exp[-2 a])
Plot[{fs[.5/x],fa[.5/x]},{x,.25,1}]

Plot[{fs[x/2]+.25/x^2,fa[x/2]+.25/x^2},{x,.25,4},PlotRange->{{0,4},{-.6,.5}}]
Plot[fs[x/2]+.25/x^2,{x,.9,4}]  

e1=Table[{1,i},{i,1,31}]
e2=Table[{1,i},{i,1,31}]
et1=Table[{1,i},{i,1,31}]
et2=Table[{1,i},{i,1,31}]
Do[
ai=.25+.025*(i-1) ;
s=s1 /. FindRoot[1/ai==s1(Tanh[s1]+1),{s1,.5/ai}] ;
es=-.5*(2*s*ai)^2 ;
e1[[i]]={ai,es}; 
et1[[i]]={1/ai,es+.25*ai^2},{i,1,31}]
ListPlot[e1]
ListPlot[et1]

Do[
ai=.25+.025*(i-1) ;
s=s1 /. FindRoot[1/ai==s1(Coth[s1]+1),{s1,.67(1/ai-1)}] ;
ea=-.5*(2*s*ai)^2 ;
e2[[i]]={ai,ea}; 
et2[[i]]={1/ai,ea+.25*ai^2},{i,1,30}]
e2[[31]]={1,0}
et2[[31]]={1,+.25*ai^2}
ListPlot[e1]
ListPlot[e2]
Show[%,%%,PlotRange->All]
Show[%,%3,PlotRange->All]

ListPlot[et1]
ListPlot[et2]
Show[%,%%,PlotRange->All]
Show[%4,%,PlotRange->{{0,4},{-.6,.5}}]

Plot[fs[.5*x]+.25/x^2,{x,.9,4}]
ListPlot[et1]
Show[%,%%,PlotRange->All]

GraphicsGrid[{{%21},{%25},{%29}}]                                      

Export["delta2_RR2.eps",%]                                              


Note parametric formula for exact E(a):
(a,E)=(s(Coth[s]+1),-(2/(Coth[s]+1))^2/2) odd
(a,E)=(s(Tanh[s]+1),-(2/(Tanh[s]+1))^2/2) even


ParametricPlot[{s(Coth[s]+1),-(2/(Coth[s]+1))^2/2},{s,.001,2}]
ParametricPlot[{s(Tanh[s]+1),-(2/(Tanh[s]+1))^2/2},{s,0,2}]

