> 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
>