How to code a formula with double integration

Hello, I've tried coding the equation in this figure using the given parameters and the functions. First i defined Fh, fU and k seperately and then called them into the integral part.And when defining them, i had to define u and epsilon as well. But i think as we have to do the double integration for u and epsilon, that two variables cannot be assigned. And that gives errors when defining Fh, fU and k. So to solve this then i tried to code the whole eqaution at once but it gave more errors.

%Parameters
n=10;
m=5;                %M Average number of users per hotspot
P_s=10^(20/10);     %BS transmit powers
P_m=10^(50/10);
alpa_L=2.0;         % Path-loss exponent       
alpa_NL=3.3;
R_s=30;             % SBS 
R=200;              % MBS
Pr=50;              %Rate threshold
beta=10^(70/10);            %Path loss at 1 m
G=10^(18/10);               %Main lobe gain
W=600e6;              % Bandwidth
NoW=10^((10*log10(W)-173+10)/10)*1000e6;    %Noise power
m_L=2;              % Parameter of Nakagami distribution
m_NL=3;
miu=170;              % LOS rage constant

Et=0:0.05:1;  
% Et=0.3% backhaul and acces distribution
W_a=(1-Et).*W;

x=(R-R_s).*rand(10,1);      %location of the node
xo=(2*pi)*rand(10,1);
u=(R_s)*rand(10,1);         %location of the user
Ep=(2*pi)*rand(10,1);

k_p=(P_s/P_m).^(1/alpa_L);
t=t_calc(R_s,(k_p*x/(1-k_p^2)),(k_p^2*x/(1-k_p^2)));
for k=0:1:5
    theta_3 = ((2.^(Pr.*(t+k))./W_a)-1)
    %theta_3=0.3;
end
Fun_k=sqrt(x.^2 + u.^2 + 2.*x.*u.*cos(Ep));

F_l=(((Fun_k.^alpa_L).*beta.*NoW.*theta_3)/(P_m.*G))*10e-72;
F_11=real(F_l);
F_h=gampdf(F_11,m_L);
F_nl=(((Fun_k.^alpa_NL).*beta.*NoW.*theta_3)/(P_m.*G))*10e-74;
F_n11=real(F_nl);
F_nh=gampdf(F_n11,m_NL);
if (0<x<R_s)
    fU=x/pi.*R_s;
else
    fU=0;
end
H=x.*k_p.*((1-(k_p.^2.*sin(Ep).^2)).^0.5+k_p.*cos(Ep))/(1-k_p.^2);
h=min(H);
u_max=min(h,R_s);
% fU= piecewise(0<x<R_s,x/pi.*R_s.^2, R_s<x,0,1);
% S=(((Pr.*Fun_k.*F_h)+((1-Pr.*Fun_k).*F_nh)).*fU./2*pi);
% P_cm=integral2(S,0,2*pi,u_max,R_s)
S=@(E,u)(((Pr.*(sqrt(x.^2 + u.^2 + 2.*x.*u.*cos(Ep))).*F_h)+((1-Pr.*(sqrt(x.^2 + u.^2 + 2.*x.*u.*cos(Ep)))).*F_nh)).*fU./2*pi);
P_cm=integral2(S,0,2*pi,0,R_s)
function T=t_calc(r1,r2,d)
T=sqrt(d+r1+r2).*sqrt(d+r1-r2).*sqrt(d-r1+r2).*sqrt(-d+r1+r2);
end

[Torsten]

And where is what you've tried so far ?

[Torsten]

E does not appear in the definition of S.

Given scalar values for E and u (which I interprete should be x and y) should result in a scalar value for S(E,u). This is not the case in your function definition: e.g. x and E_p are vectors.