Integration MATLAB QUESTION

Dear friends, i am really in trable, i am trying to find the specific attenuation due to rain but i have problem in my code and i am asking you your help. My problem is in the last command which is an integration operation, and i am sure you will do your best for helping me. The link below is for the equations and it's followed by the code.

../imgqa/eboard/Antenna/rf-utxnlysoi2s.png

radius = 1;

nMax = 40; % maximum mode number

No=8*10^3; %m^ -4

R1=140; %rainfall rate=1mm/hr

AS1=8.2*R1^-0.21;

N1D=No*exp(-AS1*radius*1e-3);

% mode numbers

mode = 1:nMax;

frequency = 6e9;

% speed of light

c = 299792458.0;

lambda = c / ( frequency ) ;

for n=1:10;

n2 = (2*n+1);

% radian frequency

w = 2.0*pi*frequency;

% wavenumber

k = w/c;

% conversion factor between cartesian and spherical Bessel/Hankel function

s = sqrt(0.5*pi/(k*radius));

% compute spherical bessel, hankel functions

[J(mode)] = besselj(mode + 1/2, k*radius); J = J*s;

[H(mode)] = besselh(mode + 1/2, 2, k*radius); H = H*s;

[J2(mode)] = besselj(mode + 1/2 - 1, k*radius); J2 = J2*s;

[H2(mode)] = besselh(mode + 1/2 - 1, 2, k*radius); H2 = H2*s;

% derivatives of spherical bessel and hankel functions % Recurrence relationship, Abramowitz and Stegun Page 361

kaJ1P(mode) = (k*radius*J2 - mode .* J );

kaH1P(mode) = (k*radius*H2 - mode .* H );

% Ruck, et. al. (3.2-1)

An = -((1i).^mode) .* ( J ./ H ) .* (2*mode + 1) ./ (mode.*(mode + 1));

% Ruck, et. al. (3.2-2), using derivatives of bessel functions

Bn = ((1i).^(mode+1)) .* (kaJ1P ./ kaH1P) .* (2*mode + 1) ./ (mode.*(mode + 1));

Qt = (lambda^2/2*pi)*sum(n2).*sum(real(An + Bn));

end

Co1=4.343*No;

================================================== ============

NOW SHOULD BE THE INTEGRATION OPERATION