Question: numerical integral of a complex expression involving BesselK function

Hi, I would like to numerically work out the profile of one expression (integral and BesselK function involved). But the computer was just jammed and came out nothing after long time waiting. I attach the worksheet below and maybe you could have a look at it or directly run it in your machine to check what is wrong, either my computer or the code itself. Thank you very much.

> restart;
> with(PDEtools); with(Units[Standard]); with(ScientificConstants);
[CanonicalCoordinates, ChangeSymmetry, CharacteristicQ,

CharacteristicQInvariants, ConservedCurrentTest, ConservedCurrents,

ConsistencyTest, D_Dx, DeterminingPDE, Eta_k, Euler, FromJet,

InfinitesimalGenerator, Infinitesimals, IntegratingFactorTest,

IntegratingFactors, InvariantSolutions, InvariantTransformation, Invariants,

Laplace, Library, PDEplot, PolynomialSolutions, ReducedForm,

SimilaritySolutions, SimilarityTransformation, SymmetrySolutions,

SymmetryTest, SymmetryTransformation, TWSolutions, ToJet, build, casesplit,

charstrip, dchange, dcoeffs, declare, diff_table, difforder, dpolyform,

dsubs, mapde, separability, splitstrip, splitsys, undeclare]
[*, +, -, /, <, <=, <>, =, Im, Re, Unit, ^, abs, add, arccos, arccosh, arccot,

arccoth, arccsc, arccsch, arcsec, arcsech, arcsin, arcsinh, arctan, arctanh,

argument, ceil, collect, combine, conjugate, convert, cos, cosh, cot, coth,

csc, csch, csgn, diff, eval, evalc, evalr, exp, expand, factor, floor, frac,

int, ln, log, log10, max, min, mul, normal, polar, root, round, sec, sech,

seq, shake, signum, simplify, sin, sinh, sqrt, surd, tan, tanh, trunc, type,

verify]
[AddConstant, AddElement, AddProperty, Constant, Element, GetConstant,

GetConstants, GetElement, GetElements, GetError, GetIsotopes, GetProperties,

GetProperty, GetUnit, GetValue, HasConstant, HasElement, HasProperty,

ModifyConstant, ModifyElement]
> Digits := 5;
5
> sigma := Units[Standard][`*`](6.5, Units[Standard][`^`](10, Units[Standard][`-`](19))); P := 50; k := Units[Standard][`*`](1.3, Units[Standard][`^`](10, Units[Standard][`-`](23))); T := 300; `#msub(mi("n"),mi("e"))` := Units[Standard][`*`](1.55, Units[Standard][`^`](10, 12)); `#msub(mi("E"),mn("0"))` := 5360; e := Units[Standard][`*`](1.6, Units[Standard][`^`](10, Units[Standard][`-`](19))); Q = Units[Standard][`*`](1000, e);
-19
6.5000 10
50
-23
1.3000 10
300
12
1.5500 10
5360
-19
1.6000 10
-16
Q = 1.6000 10
> `#msub(mi("n"),mi("n"))` := Units[Standard][`*`](P, Units[Standard][`/`](Units[Standard][`*`](k, T))); l := Units[Standard][`^`](Units[Standard][`*`](`#msub(mi("n"),mi("n"))`, sigma), Units[Standard][`-`](1));
22
1.2820 10
0.00012000
> lambda := Units[Standard][`^`](Units[Standard][`*`](Units[Standard][`*`](Units[Standard][`*`](e, `#msub(mi("E"),mn("0"))`), l), Units[Standard][`/`](Units[Standard][`*`](Units[Standard][`*`](Units[Standard][`*`](4, Pi), `#msub(mi("n"),mi("e"))`), Units[Standard][`^`](e, 2)))), Units[Standard][`/`](2));
(1/2)
/1 \
805.21 |--|
\Pi/
> X := Units[Standard][`+`](1, Units[Standard][`-`](sqrt(Units[Standard][`+`](1, Units[Standard][`*`](I, t)))));
(1/2)
1 - (1 + I t)
> Y := Units[Standard][`+`](Units[Standard][`*`](Units[Standard][`*`](Units[Standard][`*`](2, sqrt(Units[Standard][`+`](1, Units[Standard][`*`](I, t)))), Units[Standard][`/`](Units[Standard][`*`](I, t))), int(Units[Standard][`/`](Units[Standard][`^`](Units[Standard][`+`](1, Units[Standard][`*`](Units[Standard][`*`](I, t), Units[Standard][`+`](1, Units[Standard][`-`](Units[Standard][`^`](x, 2))))), 2)), x = 0 .. 1)), Units[Standard][`-`](Units[Standard][`/`](Units[Standard][`*`](Units[Standard][`*`](I, t), Units[Standard][`+`](1, Units[Standard][`*`](I, t))))));
(1/2) / (1/2) / t \\
I (1 + I t) |I (-(I - t) t) + arctanh|-----------------||
| | (1/2)||
\ \(-(I - t) t) //
- -------------------------------------------------------------------
(1/2)
t (I - t) (-(I - t) t)

I
+ -----------
t (1 + I t)
> Phi := Units[Standard][`*`](Units[Standard][`*`](Units[Standard][`*`](2, Q), Units[Standard][`/`](Units[Standard][`*`](Pi, l))), Re(int(Units[Standard][`*`](Units[Standard][`/`](Units[Standard][`+`](1, Units[Standard][`*`](Units[Standard][`^`](Units[Standard][`*`](l, Units[Standard][`/`](lambda)), 2), Y))), BesselK(0, Units[Standard][`*`](Units[Standard][`*`](r, Units[Standard][`/`](l)), sqrt(Units[Standard][`*`](Units[Standard][`+`](Units[Standard][`^`](t, 2), Units[Standard][`*`](Units[Standard][`^`](Units[Standard][`*`](l, Units[Standard][`/`](lambda)), 2), X)), Units[Standard][`/`](Units[Standard][`+`](1, Units[Standard][`*`](Units[Standard][`^`](Units[Standard][`*`](l, Units[Standard][`/`](lambda)), 2), Y)))))))), t = 0 .. infinity)));
Warning, computation interrupted
> plot(Units[Standard][`*`](Units[Standard][`*`](Units[Standard][`*`](2, Q), Units[Standard][`/`](Units[Standard][`*`](Pi, l))), Re(int(Units[Standard][`*`](Units[Standard][`/`](Units[Standard][`+`](1, Units[Standard][`*`](Units[Standard][`^`](Units[Standard][`*`](l, Units[Standard][`/`](lambda)), 2), Y))), BesselK(0, Units[Standard][`*`](Units[Standard][`*`](r, Units[Standard][`/`](l)), sqrt(Units[Standard][`*`](Units[Standard][`+`](Units[Standard][`^`](t, 2), Units[Standard][`*`](Units[Standard][`^`](Units[Standard][`*`](l, Units[Standard][`/`](lambda)), 2), X)), Units[Standard][`/`](Units[Standard][`+`](1, Units[Standard][`*`](Units[Standard][`^`](Units[Standard][`*`](l, Units[Standard][`/`](lambda)), 2), Y)))))))), t = 0 .. infinity))), r = Units[Standard][`-`](0.1e-2) .. 0.1e-2);
Warning, computation interrupted
>

Download pic_potential.txt