(* Input *)
(*kin = 20.;*)
mass= 939.;
range = 3;

sps={{Dashing[None],RGBColor[0,0,0]}, {Dashing[Large],RGBColor[1,0,0]},{Dashing[Medium],RGBColor[0,0,1]},{Dashing[None],RGBColor[0,1,0]},{Dashing[None],RGBColor[0,1,1]}};

kOut[e_] := ((Sqrt[2*mass*e])/197.327)

kIn[e_,U0_] := ((Sqrt[2 mass (e- U0)])/197.327)

ang = 0;
Plot[{e*SphericalBesselJ[ang,e],Sin[e-ang*(Pi/2)]},{e,0,50},PlotStyle->sps]

Plot[{-e*SphericalBesselY[ang,e],Cos[e-ang*(Pi/2)]},{e,0,50},PlotStyle->sps]

δ[e_,ang_,r_,u0_]:= ArcCot[(SphericalBesselY[ang,kOut[e] r]D[SphericalBesselJ[ang,kIn[e,u0]r],r] - SphericalBesselJ[ang,kIn[e,u0]r]D[SphericalBesselY[ang,kOut[e] r],r])/(SphericalBesselJ[ang,kOut[e]r]D[SphericalBesselJ[ang,kIn[e,u0]r],r] - SphericalBesselJ[ang,kIn[e,u0]r]D[SphericalBesselJ[ang,kOut[e]r],r]) ]

myDelta[e_,ang_,u0_]:= Simplify[δ[e,ang,r,u0]/.r->range]

myU = -20;
myAng = 0;
Plot[{myDelta[e,myAng,-2],myDelta[e,myAng,-5],myDelta[e,myAng,-7]},{e,0,40},PlotRange->All,PlotStyle->sps]

CS[e_,ang_,u0_]:=(1/100)((4 π )/((kOut[e])^2)) (2ang +1) (((Sin[myDelta[e,ang,u0]]))^2)

csAng = 1;
Plot[{CS[ee,csAng,-5],CS[ee,csAng,-10],CS[ee,csAng,-15],CS[ee,csAng,-20]},{ee,0,40},PlotRange->All,PlotStyle->sps]

A[e_,ang_,u0_]:=((kOut[e] SphericalBesselJ[ang,kOut[e]range]Cos[myDelta[e,ang,u0]] -kOut[e] SphericalBesselY[ang,kOut[e] range]Sin[myDelta[e,ang,u0]])/(kIn[e,u0]SphericalBesselJ[ang,kIn[e,u0]range]))

fIn[e_,r_,ang_,u0_]:=A[e,ang,u0]kIn[e,u0]r SphericalBesselJ[ang,kIn[e,u0]r]

fOut[e_,r_,ang_,u0_]:=kOut[e]r SphericalBesselJ[ang,kOut[e]r]Cos[myDelta[e,ang,u0]] - kOut[e]r SphericalBesselY[ang,kOut[e]r]Sin[myDelta[e,ang,u0]]

myE = 30;
myU = -20;
myAng = 1;
Plot[{fIn[myE,rr,myAng,myU],fOut[myE,rr,myAng,myU]},{rr,0.1,15},PlotRange->All,PlotStyle->sps]

myR = 15;
myU = -15;
myAng = 1;
Plot[{fIn[ee,myR,myAng,myU],fOut[ee,myR,myAng,myU]},{ee,0.1,40},PlotRange->All,PlotStyle->sps]

fOut1[e_,r_,ang_]:=kOut[e] r  SphericalBesselJ[ang,kOut[e]r]

TMatriX[e_,u_,ang_]:=(1/(kIn[e,u]kOut[e]))NIntegrate[ u fIn[e,rr,ang,u]fOut1[e,rr,ang] ,{rr,0,range}]

Clear[TMatriX0]
(* L=0 *)
TMatriX0[e_,u_,aa_,bb_]:=Integrate[Sin[aa rr]Sin[bb rr],rr]
Print[Simplify[TMatriX0[e,u,aa,bb]]]

tmpT0[e_,u_,aa_,bb_,rr_]:= ((u A[e,0,u])/(aa bb  ))((bb Cos[bb rr] Sin[aa rr]-aa Cos[aa rr] Sin[bb rr])/((aa^2)-(bb^2)))

u0 = -20;
T0[e_]:= tmpT0[e,u0,kIn[e,u0],kOut[e],3]-tmpT0[e,u0,kIn[e,u0],kOut[e],0]

Plot[((T0[e])^2),{e,0,50}]

(* L=1 *)

TMatriX1[e_,aa_,bb_]:=Integrate[(((Sin[aa rr])/(aa rr)) - Cos[aa rr])(((Sin[bb rr])/(bb rr)) - Cos[bb rr]),rr]

Print[Simplify[TMatriX1[e,aa,bb]]]

tmpT1[e_,u_,aa_,bb_,rr_]:= ((u A[e,0,u])/(aa bb  ))(((aa^2) bb rr Cos[bb rr] Sin[aa rr]-(aa (bb^2) rr Cos[aa rr]+((aa^2)-(bb^2)) Sin[aa rr]) Sin[bb rr])/(bb ((aa^3)-aa (bb^2)) rr))

u0 = -20;
T1[e_]:= tmpT1[e,u0,kIn[e,u0],kOut[e],3]-tmpT1[e,u0,kIn[e,u0],kOut[e],0.0001]

Plot[((T1[e])^2),{e,0,50}]