Look at K:= your integrand by plotting:

  plot3d(K, t=0 .. 70, dd=0 .. 100, axes=boxed);

That lets you guess what happens: for any dd there is just 
1 maximum below 900 and may flatten out indeed:

  diff(K, t); 
  solve(%, t);
  eval(K, t=%);
                                  17
                                 870

So the maximum is in t=17 with value = 870

Thus all integrals are less than 870 * length = 60900

Which is what you also get for dd --> +oo

But your V = 57002400 is larger.

I think there is a fast way with timing below 1/2 second (instead of
some thousands) allowing to plot the roots vs. intensity over the
range 0 ... 10^8.

I have not tried to give a compiled version as well.

For that I used Brent's root finder (having some idea about needed
'bracketting' from experiments allready done, but not needed) and
optimized evaluation of the system.

For finer results one probably needs better data (for example ep0,
is it 1/(4*3.14) or is it 1/(4*Pi) ?), and higher precision.


Find my links appended (the board does not upload my files properly).


http://axelvogt.de/temp/MP_fsolve_2dim_sytem_dilogs_fast.mws
http://axelvogt.de/temp/MP_fsolve_2dim_sytem_dilogs_fast.pdf

The 2-dim system are all the same and only differ by the 'intensity',
where the 2nd equation is the same in any case (an originator should
state that in a clear way).

http://www.mapleprimes.com/questions/201767-Fsolve-Cannot-Solve-This-Equations-And
http://www.mapleprimes.com/questions/201804-Reliability-Result-Of-DirectSearch
http://www.mapleprimes.com/questions/201813-DirectSearch-Result

I set up the system using the symbolics only and keep data separately.


For 14250 <= intensity one can estimate the solution by the following
up to an relative error of 2*e-9, with almost immediate result:

  tmp3 := (-aa*enph^3-6*u)*exp(-lambda-c)-exp(-c-b)*aa*enph^3+
    (-aa*enph^3+6*u)*exp(lambda-b)-aa*enph^3;
  approx_eq2 := c = -1/2*(2*((-Pi^2+b^2-2*impurity/cc)^2+
    4*Pi^2*b^2)^(1/2)-2*Pi^2+2*b^2-4*impurity/cc)^(1/2);

  'eval(tmp3, approx_eq2)';
  eval[recurse](%, allData);         # all the actual data 
  fsolve(%,b);

If one uses the more expensive

  eval[recurse](tmp3, allData);      # all the actual data:
  phi:=unapply(%, b,c);
  initialGuess:=fsolve(phi(b,'C'(b)), b=bMin .. infinity);

the result also comes up quite quick, error ~ 2*e-12

Where the numerical solution for eq2 can be written as

  C:= b -> 
    fsolve(polylog(2,-exp(-c)) = polylog(2,-exp(b))+1/25*Pi, c=-b);

and bMin:= -2.04299751639794 (solving polylog(2,-exp(b))+1/25*Pi).

PS: I used lambda instead of l which does not look like 1. Digits:=15.

Besides Markiyan's reply here is some intrinsic, but ad hoc way.

Look at sys = {F,G}, where G ~ -polylog(2,-exp(b))+polylog(2,-exp(-c))- Pi/25
has a zero locus = curve "like" a (partially) defined function ('implicitplot'
shows that). 

Use a numerical procedure C to define that function. Since 'implicitplot' also
indicates some linear behaviour use that to provide a good initial guess (for
speed up). I did that by a series approach.

For the system now one has to solve F(b, C(b)). Maple can do it (you mean Reals?).

http://www.axelvogt.de/temp/MP_fsolve_2dim_sytem_dilogs.mws
http://www.axelvogt.de/temp/MP_fsolve_2dim_sytem_dilogs.pdf
 

Here is a somewhat brute suggestion for your second question, simple case

select_order:=proc(DE, k::nonnegint)
local P, Delta;
  P:=convert(DE, diff);
  P:=eval(P, diff = '(f,x) ->Diff(f,x)*Delta' );
  P:=eval(P, Diff=diff);
  coeff(P, Delta, k);
end proc; 

Tested it (DE should be free of Delta, or refine the code (*) ):

  tst:= eq;
  [seq(select_order(%, n), n=0 .. 5)];
 
  tst=convert(%, `+`): # recollect
  is(%);

which gives the desired and recollection works (giving 'true').

It works for more complicated expressions, try it for
 
  tst:=(1+eq)*x;
  tst:=diff((1+eq)*x, x): tst:=expand(tst):

PS & edited: something like DE <---> Polynom probably is covered
by packages on differential algebra (or physics?) in a sound way.

(*) I now made Delta local, as suggested by Carl

evalf[50](10.2^20);

It first evaluates the inner term. It is a matter of taste whether the next command should append zeros (since we are in decimal representation)

convert(%, radical) will do it

And cos(...) is accurate, of course, you want algebraic expressions for that.
But not all users want to switch from trigonometrics to radicals automatically.

PS: please post code, not pictures

H1:=proc(y, eps)                # similar: H = without an 'epsilon' 
  evalf( Int('x -> (-12*y^2+1)/(4*y^2+1)^3 * ln(abs(Zeta(x+y*I))) ',
    1/2 .. infinity, method = _d01amc, epsilon=eps) );
end proc;

`task` = Int( 'H'(y), y=0 .. infinity);

Plotting near zero looks nice, but it has a wild behaviour and does
not even exhibit a good decay, try on y = 50 .. 100 as example.

Now integrate on y = 2^k .. 2^(k+1) and with reasonable error control
I final get the result as 0.2378564, 6 correct decimals (without any
proof, why it is justified to stop that sum ...)

Find it appended (if Maple's upload works).

MP_integral_2dim_zeta.mws

MP_integral_2dim_zeta.pdf

Für welche Variablen? Es sind wenigstens 4: mu,nu, lambda, alpha - welche ist die parametrische Variable?

Download plot-it.mw

After re-writing it is a polynomial in y of degree = 3, Maple can solve that and the first of the three solutions is the desired real function. If you want the complex solutions as well then you may use abs(...).

Over the complex numbers the situation is similar: only along the axes,
in boxes like the real case as the following plot suggests:

m,n:= 3,5;
myRange:= x= -m*Pi .. m*Pi, y= -n*Pi .. n*Pi;
tan(x+y*I)-tanh(x+y*I);
plot3d(abs(%), myRange, view=-1 .. 1, axes=boxed, orientation=[90,0]);

You may use tan(z)+tan(z*I)*I = tan(z)-tanh(z) to check symmetries.

Initial guess for extreme situations is a task to be done ...

If I understand correctly than the following is assumed:

  assume(0<t, 1<B, 0<eta, 0<h, 0<r, 0<H); # I will use r and H later

You can rewrite your integral and do it in steps

  f:=B*1/t*(t/eta)^B*exp(-(t/eta)^B);
  'Int(f*(t-H), t)'; # = yours, indefinite, setting H = n*h

  'Int(f*H, t)';
  FH:=value(%);

                                      / t \B
                        FH := -H exp(-|---| )
                                      \eta/

Note that this is the same as changing variables and changing back:

  'Int(f*H, t)';
  IntegrationTools:-Change(%, t/eta=r, r);
  value(%);
  eval(%, r=t/eta); # change back

Now do the other

  'Int(f*t, t)';
  Change(%, t/eta=r, r); #expand(%);
  value(%):
  
  # write it as more common functions
  convert(%, hypergeom): convert(%, StandardFunctions):  
  
  # beautify that
  simplify(%): simplify(%, size):
  collect(%, [GAMMA, exp]):
  
  Ft:=eval(%, r=t/eta); # change back

I omit the lengthy output, it is in terms of Gamma function and exp
(and it can be simplied a bit more).


Now your anti-derivative is Ft - eval(FH, H=n*h) from which you can
go on, hoping the fundamental theorem applies (continuity needed).


For the objective function (as written down by acer) you 'only' need
to devide by eval(FH, H=1) to go on.

Hope that helps, an upload is not possible for me here (shameshameshame)

I think that MMA is not correct:

Changing a*r = x shows that 'a' can be ignored. The order can be reduced:

  Int(x*BesselI(1,x)*BesselK(1,x), x);
  IntegrationTools:-Parts(%,x*BesselI(1,x)): expand(%);

      -BesselK(0, x) x BesselI(1, x)

             /
            |
         +  |  BesselK(0, x) BesselI(0, x) x dx
            |
           /

For the latter Wolframalpha gives an answer

  Mma:= 
  convert( "MeijerG[{{3/2}, {}}, {{1, 1}, {0}}, x, 1/2]/(4 Sqrt[Pi])", FromMma);

which I did not trust, but one can correct it:

  convert(Mma, `0F1`):
  eval(%, x=x^2): combine(%) assuming 0<x:
  
  Diff(%,x):
  expand(value(%)) = %: simplify(%, size);

                                  d        2
  BesselK(0, x) BesselI(0, x) x = -- (1/2 x
                                  dx

        (BesselI(0, x) BesselK(0, x) + BesselK(1, x) BesselI(1, x)))

Now use that to get it for order 1:

  -BesselK(0,x)*x*BesselI(1,x) + 
    1/2*x^2*(BesselK(1,x)*BesselI(1,x)+BesselI(0,x)*BesselK(0,x)):
 
  Diff(%,x):
  expand(value(%)) = %: simplify(%, size);

                                  d
  x BesselI(1, x) BesselK(1, x) = -- (1/2 (
                                  dx

        (x BesselI(0, x) - 2 BesselI(1, x)) BesselK(0, x)

         + x BesselI(1, x) BesselK(1, x)) x)

Now plotting that solution is quite different from the MeijerG which
was initially asked (and not only by some constant).

PS: the integral is numerically quite easy (as the poster asked for
Matlab) and perhaps even more simply than a symbolic version.

No sheet attached, since the upolader does not work (neither in FF nor IE)

Edited: here are the files through external links
http://www.axelvogt.de/temp/MP_integral_BesselI1K1.mws
http://www.axelvogt.de/temp/MP_integral_BesselI1K1.pdf

I think there is a formal result for that ugly task and I sketch it.

  Int(cosh(t)/(cosh((17/15)*t)+cosh(t)), t = 0 .. infinity);
  Change(%, t=-I*z*15, z); expand(%): normal(%);
  subs(infinity=oo, %);
  Change(%, z=I*t, t);
  combine(%); #expand(%);
  #simplify(%);
  subs(oo=infinity, %): sort(%);
  
  Change(%, 16*t = Pi*x, x); expand(%): sort(%);
  
  GetIntegrand(%);
  
  A,B:=numer(%), denom(%);
  LA:=[op(A)];

That shows you are after Int( cosh / cosh(x) ), I lost a 'minus':


  [Int(-2/cosh(Pi*x)*cosh(1/8*Pi*x), x = 0 .. infinity), 
   Int( 2/cosh(Pi*x)*cosh(1/4*Pi*x), x = 0 .. infinity),
   Int(-2/cosh(Pi*x)*cosh(3/8*Pi*x), x = 0 .. infinity), 
   Int( 2/cosh(Pi*x)*cosh(1/2*Pi*x), x = 0 .. infinity),
   Int(-2/cosh(Pi*x)*cosh(5/8*Pi*x), x = 0 .. infinity),
   Int( 2/cosh(Pi*x)*cosh(3/4*Pi*x), x = 0 .. infinity),
   Int(-2/cosh(Pi*x)*cosh(7/8*Pi*x), x = 0 .. infinity),
   Int( 1/cosh(Pi*x),                x = 0 .. infinity)];

   evalf(%); convert(%, `+`); 15/32*Pi*%: evalf(%);

                          -5.21062483300027


I vaguely remembered (from the usenet)some 'stuff' of Ramanujan, his
task was cos(Pi*r*x)/cosh(Pi*x)*exp(-Pi*w*x^2) as integrand and for
r=I*something, w=0 this is your task. Filed papers cover 0 < w, so
they do not apply. However they say 'fourier is the plan to attack'
and that works for your special case.

  with(inttrans):

Then observe 2 things

  2*sqrt(Pi/2)*fouriercos( f(x), x , s):
  %=eval(%, [f = 'xi -> 1/cosh(Pi*xi)']);

             1/2   1/2                              1
            2    Pi    fouriercos(f(x), x, s) = ---------
                                                cosh(s/2)

  2*sqrt(Pi/2)*fouriercos( f(x), x , s);
  convert(%, Int):
  %%=eval(%, [f = 'xi -> 1/cosh(Pi*xi)', s=I*k/8]);
 
                                            infinity      k x
                                           /         cosh(---)
    1/2   1/2                             |                8
   2    Pi    fouriercos(f(x), x, s) = 2  |          ---------- dx
                                          |          cosh(Pi x)
                                         /
                                           0

Using all for the above list should give a formal result. Sorry
for being too lazy to write down the 8 terms.

PS: being non lazy it is
result := -15/32*Pi*
  (1/2-1/cos(1/16*Pi)+1/cos(1/8*Pi)-1/cos(3/16*Pi) +
   2^(1/2)-1/cos(5/16*Pi)+1/cos(3/8*Pi)-1/cos(7/16*Pi))

 restart;
 f:= phi -> exp( - (phi-phi0)^2/2/sigma^2);

                                                  2
                                      (phi - phi0)
                 f := phi -> exp(-1/2 -------------)
                                              2
                                         sigma

 inttrans[fourier](f(phi), phi, omega); #%  assuming 0 < sigma, phi0::real;

              2                       2
          phi0                     phi      phi0 phi
   exp(- --------) fourier(exp(- -------- + --------), phi, omega)
                2                       2         2
         2 sigma                 2 sigma     sigma

 convert(%, Int):
 value(%);

             2    /{
         phi0     |{
  exp(- --------) |{
               2  |{
        2 sigma   \{

                        2          2
            (omega sigma  I - phi0)    1/2        1            1/2
        exp(------------------------) 2    csgn(-----) sigma Pi    ,
                           2                    sigma
                    2 sigma

               1
        csgn(------) = 1
                  2
             sigma

                            \
                            |
        infinity , otherwise|
                            |
                            /

 simplify(%) assuming 0 < sigma; # , phi0::real;
 expand(%): combine(%, exp);

                                 2               1/2         1/2
     exp(1/2 I omega (omega sigma  I - 2 phi0)) 2    sigma Pi


                                     2      2
          1/2         1/2       omega  sigma
         2    sigma Pi    exp(- ------------- - omega phi0 I)
                                      2