how to simulate a function

i'm trying to simulate this function "pd" against "lnd" by accept reject monte carlo method but it takes forever and nothing comes out

M=8;
N=32;
K=14;
c=2;
s=1;
T=18.623390280430364437245197384232;
lnd=0:5:30;
pd1=(2*K*nchoosek(N/2,K)*((gamma(K+i)*gamma(N-K+i+1+T./(1+lnd)))./(gamma(N+1+T./(1+lnd)))-...
  (symsum(nchoosek(N/2,i)*(gamma(K+i)*gamma(N-K-i+1+T./(1+lnd)))./(gamma(N+1+T./(1+lnd))), i, K, (N/2)))));
%lnd=0:5:30;
pd=(symsum(nchoosek(M,i)*(pd1).^(i).* (1-(pd1)).^(M-i), i, s, M));
O=10^7;
Mx=1.2*max(pd);
lndi=10*rand(1,O);
pj=10*rand(1,O);
pd1i=2*K*nchoosek(N/2,K)*((gamma(K+i)*gamma(N-K+i+1+T./(1+lndi)))./(gamma(N+1+T./(1+lndi)))-...
    (symsum(nchoosek(N/2,i)*(gamma(K+i)*gamma(N-K-i+1+T./(1+lndi)))./(gamma(N+1+T./(1+lndi))), i, K, (N/2))));
f_lndi=(symsum(nchoosek(M,i)*(pd1i).^(i)* (1-(pd1i)).^(M-i), i, s, M));
comp=pj*Mx;
result=f_lndi>=comp;
index=find(result);
lnd_new=lndi(1,index);
[pd1,lnd1]=hist(lnd_new,100);
dlnd=(lnd1(2)-lnd1(1));
h=(pd1/length(lnd_new))/dlnd;
plot(lnd1,h,'b.');

Answer 1

[Walter Roberson]

As I told you in https://www.mathworks.com/matlabcentral/answers/342612-why-does-it-take-forever-to-solve-a-function, do not use symsum() for adding definite symbolic terms.