function [dn]=numdist(z,n)

% k and l refer to the array indexes in the D and L
% variables, specifying which of the DL combinations
% we are solving for on this function call

global d0 mfr T Po MW mu smf rhos kb R rhog Vo SECT rhomix nu startnum D L k l m
global uz A funcall

funcall=funcall+1;
dn=zeros(SECT,1);
sums=zeros(3,1);


    disp(['Function call ', num2str(funcall),' z value ', num2str(z)...
        ,' of ',num2str(L(m)),'.'])

for i=1:SECT   
   if (i-2 >= 1)
        for j=1:i-2
          sums(1)=sums(1)+beta(i-1,j)/2^(i-j-1)*n(j);
        end
   end
   if (i-1 >= 1)
         for j=1:i-1
           sums(2)=sums(2)+beta(i,j)/2^(i-j)*n(j);
         end
         dn(i)=dn(i)+ n(i-1)*sum(1)+ 0.5*beta(i-1,i-1)*n(i-1)^2-n(i)*sums(2)-gamma(i)*n(i);;
   end 
   if (i+1 <=SECT)
       	dn(i)=dn(i)+2*gamma(i+1)*n(i+1);
        for j=1:SECT
             sums(3)=sums(3)+beta(i,j)*n(j);
        end
        dn(i)=dn(i)-n(i)*sums(3);
   end   
   dn(i)=dn(i)/uz(k);
   sums=zeros(3,1);
end