Mp=938.272; (*mass of proton, MeV*)
Mn=939.565;(*mass of neutron, MeV*)
Ebind1=3.925 ;(*binding energy ((of )^8)He, MeV*)
Esep1=24.81432; (*1p separation energy ((for )^8)He, MeV*)
p=197.327;
mass1=821;
range1=3.735;

fIn[q_,r_,ang_]:=q r SphericalBesselJ[ang,q r]
fOut[k_,r_,ang_]:=((((2/Pi)k r))^(1/2)) BesselK[ang+(1/2),k r]
dfInR[q_,r_,ang_]:=D[fIn[q,r,ang],r]
dfIn[q_,ang_]:=Simplify[dfInR[q,r,ang]/.r->range1]
dfOutR[k_,r_,ang_]:=D[((((2/Pi)k r))^(1/2)) BesselK[ang+(1/2),k r],r]
dfOut[q_,ang_]:=Simplify[dfOutR[q,r,ang]/.r->range1]

Plot[((fIn[((Sqrt[2 mass1 (-Esep1+U)])/p),range1,0])/(dfIn[((Sqrt[2 mass1 (-Esep1+U)])/p),0]))-((fOut[((Sqrt[2 mass1 Esep1])/p),range1,0])/(dfOut[((Sqrt[2 mass1 Esep1])/p),0])),{U,0,50}]

myU1=U/.FindRoot[((fIn[((Sqrt[2 mass1 (-Esep1+U)])/p),range1,0])/(dfIn[((Sqrt[2 mass1 (-Esep1+U)])/p),0]))-((fOut[((Sqrt[2 mass1 Esep1])/p),range1,0])/(dfOut[((Sqrt[2 mass1 Esep1])/p),0]))==0,{U,60}]

myCoeff1=coeff/.FindRoot[fIn[((Sqrt[2 mass1 (-Esep1+myU1)])/p),range1,0]-coeff *fOut[((Sqrt[2 mass1 Esep1])/p),range1,0]==0,{coeff,1}]

myNorm1=A/.FindRoot[(A^2) Integrate[((fIn[((Sqrt[2 mass1 (-Esep1+myU1)])/p),r,0])^2), {r,0,range1}]+(((myCoeff1*A))^2) Integrate[((fOut[((Sqrt[2 mass1 Esep1])/p),r,0])^2), {r,range1,Infinity}]==1, {A,0.2}]

Plot[Piecewise[{{myNorm1 fIn[((Sqrt[2 mass1 (-Esep1+myU1)])/p),r,0],r<range1},{myNorm1 myCoeff1 fOut[((Sqrt[2 mass1 Esep1])/p),r,0],r>range1}}],{r,0,10}]

(*R для гелия8 == 1.688 фм *)

(2/3)Sqrt[Integrate[(myNorm1^2)((fIn[((Sqrt[2 mass1 (-Esep1+myU1)])/p),r,0])^2)(r^2),{r,0,range1}]+Integrate[(myNorm1^2) (myCoeff1^2) ((fOut[((Sqrt[2 mass1 Esep1])/p),r,0])^2)(r^2),{r,range1,Infinity}]]

Integrate[Sin[a x] Exp[-b x],{x,r,Infinity},Assumptions->a>0 && b>0]

Integrate[Sin[a x] Sin[b x],x,Assumptions->a>0 && b>0]

PsiP[q_]:= (1/q)(Sqrt[(1/(p ))])myNorm1 (((Sin[((q/p)-((Sqrt[2 mass1 (-Esep1+myU1)])/p)) range1])/(2 ((q/p)-((Sqrt[2 mass1 (-Esep1+myU1)])/p))))-((Sin[((q/p)+((Sqrt[2 mass1 (-Esep1+myU1)])/p)) range1])/(2 ((q/p)+((Sqrt[2 mass1 (-Esep1+myU1)])/p))))+myCoeff1(((E^(-((Sqrt[2 mass1 Esep1])/p) range1)) ((q/p) Cos[((q range1)/p)]+((Sqrt[2 mass1 Esep1])/p) Sin[((q range1)/p)]))/(((q^2)/(p^2))+(((((Sqrt[2 mass1 Esep1])/p)))^2))))

(* Проверим нормировку должно быть (Pi/2)*)
NIntegrate[(((q PsiP[q]))^2),{q,0,Infinity}]

Plot[PsiP[q],{q,0,1000},PlotRange->All]