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

You can try the following (the integrand is rational in C):

EQ1 := convert(eq1, rational); # to avoid floating point numbers

lhs(EQ1); 
V1:=value(%) assuming 0<C;
V1:=convert(V1, piecewise, C);

After some time you get an answer. Plot it (in C = 0 ... 1) to see if
you have any chance. You do have a good chance. Howver 'Solve fails
(a bug?). But fsolve will do.

Dito for V2 with C < 0 (I have not checked symmetry).

For V1 or V2=-1140 I get C = +- 0.100907264887074

I think the limit is 1 by the initial value theorem for Laplace transforms.

The 2nd integral works, while the 1st does not. But you forgot to define parameters which is

not n = 1 but n:= 1 (with double point, and so on)

PS: I do not correct it, since handling sheets in document mode is painful

After all the discusion here is the answer (applying Gradshteyn & Ryzhik (2005), Formula 6.512.1 and the Errata for it, giving a a scaled, normalized 2F1, but use the reciprocal Gamma instead of dividing by Gamma):

it is zero, as Markiyan asserted

  Phi:= z -> piecewise(z::nonposint, Limit(1/GAMMA(zeta), zeta=z), 1/GAMMA(z)):

  setting:=mu=1, a=1500,nu=2,b=r; # for Errata 
  # versus mu=2, a=r,nu=1,b=1500 for the Book

  # show that the settings give the task
  Int(BesselJ(mu,a*x)*BesselJ(nu,b*x),x = 0 .. infinity):
  % = subs(setting, %);

  # now apply the formula in G&R
  Int(BesselJ(mu,a*x)*BesselJ(nu,b*x),x = 0 .. infinity) = 
    b^nu*a^(-nu-1)* GAMMA((mu+nu+1)/2)/GAMMA(nu+1)/GAMMA((mu-nu+1)/2)*
    hypergeom([(mu+nu+1)/2, (mu-nu+1)/2 ],[nu+1],b^2/a^2);

  ``=b^nu*a^(-nu-1)* GAMMA((mu+nu+1)/2)*'Phi'(nu+1)*'Phi'((mu-nu+1)/2)*
    hypergeom([(mu+nu+1)/2, (mu-nu+1)/2 ],[nu+1],b^2/a^2);

  subs(setting, %); eval(%); 
  value(%);

                                  = 0

  # and be sure, that the conditions the formula (Errata) are satisfied
  [b < a, 0< a, 0 < b, -1 < Re(mu+nu)]; # Errata: b < a, not a < b !
  subs(setting, %);
  map(is, %) assuming 0 < r, r < 40;

               [b < a, 0 < a, 0 < b, -1 < Re(mu + nu)]
                       [true, true, true, true]

Download MP_task_BesselJ.mws

Using Maple 17 I would do:

Sum(1/(k^2-p),k=0..infinity); value(%); discont(%, p);

                                    1/2       1/2
                         Pi cot(Pi p   ) p + p
                    -1/2 ------------------------
                                    3/2
                                   p

                                     2
                              {0, _Z1 }

Which says: as far as p is not the square of an integer
Does not look nice or easy to understand. But does it.

gamma is a numerical constant in Maple, you can not assume on it

[solve(eq,beta)];
select('q -> is(0<=q)', %) assuming 1<g;

                               2     1/2
                             (g  - 1)
                            [-----------]
                                  g

Where I would prefer words instead: 1 < q, so sqrt is real & positive
and deviding by a positive does not change it. I see no need to do
work with software here (just my opinion)

I guess: k/6/Pi, k = 0 .. 6 and period Pi? But what is the reference for the task?

Edited: I used "convert(%, sin): combine(%): numer(%): collect(%, sin);" for my plot, then it works without increasing numpoints by plot(%, x=-0.1 .. Pi+0.03,  -1 .. 1);

x^t/y^t = (x/y)^t; 

                              t
                             x          t
                            ---- = (x/y)
                              t
                             y

subs(x=1, y=-1, t=1/2, %);

                                -I = I

Here is a different approach, working for larger numbers (but I do not know
to generalize the constraint). If in the following one sets theNumber:=100
then one gets Kitonum's answer (and the constraint is satisfied). But I do
not know, whether the solution is minimal. However it works and does not
suffer from combinatorial explosions. I show it for 200! instead of 100!.

restart: interface(version);
theNumber:=200;
theFactorial:='theNumber!'; `~~`=evalf(%);
                                             375
                     ~ = 0.788657867364791 10


ifactors(theFactorial)[2]: 
H:=sort(%);  # prime number decomposition
aF:=nops(%); # Anzahl Faktoren
                               aF := 46

F:=[seq( op(1,H[i]), i = 1 .. aF)]; # Prim-Faktoren
M:=[seq( op(2,H[i]), i = 1 .. aF)]; # Multiplizitäten
K:=[seq( F[i]^M[i],  i = 1 .. aF)]; # primary Komponenten
#'theFactorial=product(K[i], i = 1 .. aF)'; is(%);

r:='r':
sys:=NULL:
for i from 2 to aF - 1 do
  sys:= sys, (K[i]-1) + r[i]*K[i]=(K[i+1]-1) + r[i+1]*K[i+1];
end do:
 
i:=1:
sys:= sys, r[i]*K[i]=(K[i+1]-1) + r[i+1]*K[i+1]:
sys:={sys}: i:='i':

sol:='isolve(sys, t)'; #: expand(%);
sol:=%:
                        sol := isolve(sys, t)

eval(sys, sol): map(is, %);
                                {true}

eval(sys, sol):
X:=lhs( %[1]);  # that is linear in t

0 < X: solve(%): t in evalf(%);

         t  in  RealRange(Open(-0.707274961818205), infinity)

z:=eval(X, t=0); # = coeftayl(X, t=0, 0)

  z := 55779796302805917999479190314474522489737682702230310402429\
        7818973188035586150018567544726120267456829264718819404726\
        7416941228242282110716071584580903149517863596346964322709\
        3217976455085566203844368409127431658463316731949007653535\
        0330428665588959653463722356562158478360867394606691475414\
        2768977006039823370360752097105410374234232308961090041066\
        19977392256259918212890624

`~~`=evalf(z);
                                             375
                     ~ = 0.557797963028059 10

#z/(99!); evalf(%);
evalf(z / (theFactorial/theNumber));

                           141.454992363641

Now verify it is a solution

modp(z*(z+1), theFactorial);
                                  0

For a more detailled look one can use the following code(s):

#for i from 1 to aF do
#  modp(z, K[i])* modp(z+1, K[i]), 'component'[i] = H[i];
#end do; i:='i':

#for i from 1 to aF do
#  modp(z, K[i]), modp(z+1, K[i]), 'component'[i] = H[i]; #, 'elements' = K[i];
#end do; i:='i':

The last code was the reason for the above - feeding it with Kitonum's
solution z and for theFactorial = 100! and experimenting & pondering.

I always hated such notations in Math, but besides Carl's answer:

limit(delta(x), delta = 1);
                                  x

diff(delta(x), x);
                                delta

Better multiply by cosh(C) to avoid troubles.

Then note that Maple wants to find a symbolic solution (over the complex numbers) and does not find it (even using 'allvalues' will not help).

So if you want real solution you better use fsolve.

If you plot that task you will see, that there are 2 solutions:

cosh(C+cosh(C)) - 2*cosh(C);
plot(%, C= -2.5 .. 0.5);

fsolve(f, C=-2.5);
                          -2.38320860635885
fsolve(f, C= 0.5);
                          0.323074102021700

convert(%, Diff) will show what you want to see (UpperCase D)