How to numerically integrate with variable bounds

In short, I am trying to integrate:

where

So my struggle is trying to integrate with respect to g.

my code thus far is:

    gmin = 2;
    gmax = 6;
    p = 2.5;
      
    critmax = gmax.^2;
    v= 0:.05: 2.*critmax;
    syms g;
     
    in = @(v,g)  F(v,g.^2).*g.^-p;
    pow= @(v,min,max) vpaintegral(in,g,gmin,gmax);   
    

and F is defined here:

    inner = @(z) besselk(5/3,z);
    f= @(x,y) x.* integral(inner,x./y,Inf);
    F= @(x,y) arrayfun(@(x)f(x,y),x);

The error I recieve is that integral limits must be scalar or double

I know this is because I am using a symbolic g to calculate F, but I can't think of any other way to calculate this nested integral. Perhaps some double integral function I am unaware of, or a different way to represent g?

Answer 1

[Torsten]

gmin = 2;
gmax = 6;
p = 2.5;
inner = @(x,g) x*integral(@(z)besselk(5/3,z),x/g^2,Inf);
fun = @(x)integral(@(g)g^(-p)*inner(x,g),gmin,gmax,'ArrayValued',true);
x = 0.1:0.1:20;
y = arrayfun(fun,x)
plot(x,y)

[Nicolas Lamb]

Thank you, makes alot more sense