Integrate $\int_{0}^{\infty }e^{-\lambda\pi\big( k_1(r) + k_2(r) -k_3(r)\big)}f(r)dr$

Question

This is my Matlab code for the above function. Am I missing something?

big_X0=10^(-27.96/10);
P_tx_h = big_X0/(5^(-2));
alpha = 4;
big_X0=10^(-27.96/10);
delta = 2/alpha;
P_c = 1;% 
P_d = 1;
%T = 0:2:20;
lambda_c = 0.01;
lambda_d = 0.01;
mu=0.7854;
Tplot= 10.^(0/10);
f=zeros(1,11,'sym');
k1 = zeros(1,11,'sym');
k2 = zeros(1,11,'sym');
k3 = zeros(1,11,'sym');
p_c_uic = zeros(1,11,'sym');
fun = zeros(1,11,'sym');
f = zeros(1,11,'sym');
big_X= zeros(1,11);
r_dist = zeros(1,11);
%big_X=10.^((0:2:20)/10);
T = zeros(1,10);
for i = 1:10
big_X(i)=  10.^(T(i)/10);
r_dist(i) = big_X(i)/P_tx_h;
v = (Tplot.^delta).*gamma(1-delta).*gamma(1+delta);
lambda_c_eq_intf = lambda_c + (lambda_d.*(P_d/P_c).^delta).*(v./mu);
lambda_d_eq_intf = lambda_d + (lambda_c.*(P_c/P_d).^delta);  
f(i)= 2.*pi.*lambda_c.*r_dist(i).*exp(-pi*lambda_c.*r_dist(i).^2);
T(i) = T(i)+2;
end
for r = 1: 10
syms n
k1(r)= @(r_dist) (delta/(1-delta)).*Tplot.*r_dist(r).^2.* hypergeom([1, 1-delta], 2-delta,-Tplot);
k2(r)=  @(r_dist)symsum((((-1).^(n+1).*((big_X(r)./P_c).^n).*(r_dist(r).^(2+(n.*alpha))).* hypergeom([1, -delta-n],-delta-n + 1,-Tplot))/factorial(n).*((1+(n*alpha/2)))), n, 0, Inf);
k3(r)=  @(r_dist)delta.*((P_c./big_X(r)).^delta).*gammainc(delta, (r_dist(r).^(alpha).*(big_X(r)./P_c)));
fun(r) = @(r_dist) exp((-lambda_c_eq_intf.*pi).*(k1(r) + k2(r) -k3(r))).*f(r_dist);
p_c_uic(r) = integral(fun(r),0,Inf); % This is where I get the error
plot(big_X,p_c_uic)

end

@Edited the code based on fedja's comment