Try this
x=[0:.1:11]
y=[0:.1:11]
[X,Y]=meshgrid(x,y)
q=[1:1:11]
before integrating let's have a look to the function for k=1 q(k)=1
R=(X^2+Y^2).^.5
F2=(R./(1+R)).^.5
mind the points to make sure element-to-element don't go array and the result is Complex
before integrating it's also recommended to have a look to the function just in case there are negative values that may reduce the result inadvertently
F=q(k)*(X.^2+Y.^2)^.5./(1+q(k)*sqrt(X.^2+Y.^2)^.5)^.5
surf(X,Y,F2)
Now integrate, with the nested trapz you asked for:
I=trapz(y,trapz(x,F,2))
I = 1.201913622981376e+02
now you may want q to sweep:
F3=(q(k)*R./(1+q(k)*R)).^5
I3=zeros(1,length(q))
for k=1:1:length(q)
I3(k)=trapz(y,trapz(x,F3,2))
end
for the given values of q the integral is always the same
stem(I3)
var(I3)
ans = 2.221432309102369e-28
If interested on the primitive, MuPAD (Symbolic Toolbox required) shows the following integration primitive:
that as complicated as it seems, it is not. The imaginary part is null atan(i)=0 and the real part is independent of r, yielding an apparently absurd 1/0 that is the way MuPAD has to say 'don't know this way, try with limits'. But we have already seen that the function is constant for large R.
If you find this answer useful please click on the thumbs up and flag above thanks in advance John