Symbolic derivative of an arrayfun

Hello everyone,

I am trying to compute the derivative of the piecewise function F_Ma, which I called f_Ma.

infty = 1e9;   zero = 1e-9 ;  Pi = 3.1415 ;

%% PARAMETERS
type =1;
p = [1.585 1.585 10] ; 
H=100 ;  NA = 10 ;
urbLev = 3 ; 
alpha = [2, 3, 3]  ;  
tau = 0.31623 ;   sigma_n2 = 1e-12 ; 
m = [2 1 1] ;    
Rd  = 1e3 ;  
Ru = [0 0.3 0.7]' * Rd  ; 
eta = [0.9772  0.007943  1 ;  0.7943  0.01  1 ;  0.6918  0.005  1 ;  0.5888  0.0004  1]  ;   eta(:,3) = eta(:,1) ;
Sa = [4.88 9.6117 12.0810 27.304]' ;   Sb = [0.429 0.1581 0.1139 0.0797]' ; 

nsteps = 1 ;   
Z = linspace(zero, max(Ru)+Rd, nsteps) ;

%% CDF ANALYSIS
for u = 1:length(Ru)
IND{u,1} = @(z) z<= Rd-Ru(u) ;
IND{u,2} = @(z) (Rd-Ru(u) < z) & (z < Rd+Ru(u)) ;
IND{u,3} = @(z) z >= Rd+Ru(u) ;
    Betai{u} = @(z) real( (z>0).*acos( (z.^2 + Ru(u).^2 - Rd.^2) ./ ( (z>0).*2.*z.*Ru(u) + (z==0)) ) )  ;
    dBetai_dz{u} = @(z) (z>0 ).*(Ru(u).^2 - z.^2 - Rd.^2) ./ ...
                            ( (z>0) .* (z .* sqrt(4*z.^2 .* Ru(u).^2 - (z.^2 + Ru(u).^2 - Rd.^2).^2)) + (z==0) )   ;
    Betai{u} = @(z) real( acos( (z.^2 + Ru(u).^2 - Rd.^2) ./ ( 2.*z.*Ru(u) ) ) )  ;
    dBetai_dz{u} = @(z) (Ru(u).^2 - z.^2 - Rd.^2) ./ ...
                            ( z .* sqrt(4*z.^2 .*Ru(u).^2 - (z.^2 + Ru(u).^2 - Rd.^2).^2) )   ;
    zX{u} = @(beta) Ru(u).*cos(beta) + sqrt(Rd.^2 - Ru(u).^2.*sin(beta).^2) ; 
    Sigma{1} = @(z) pi*z.^2 ;
    Sigma{2} = @(z) Betai{u}(z).*z.^2 + 1/2*integral(@(beta) zX{u}(beta).^2, Betai{u}(z),2*pi-Betai{u}(z),'ArrayValued',1) ;
    Sigma{3} = @(z) pi*Rd^2 ;
dSigma{1} = @(z) 2*pi*z ; 
dSigma{2} = @(z) 2*z.*Betai{u}(z) + z.^2 .* dBetai_dz{u}(z) - zX{u}(Betai{u}(z)).^2 .* dBetai_dz{u}(z) ; 
dSigma{3} = @(z) 0 ; 
P_Ma{1} = @(z) 1./(1+Sa(urbLev)*exp(-Sb(urbLev)*(180/pi*atan(H./z)-Sa(urbLev))))  ;
P_Ma{2} = @(z) 1-P_Ma{1}(z)  ;
for M = 1:2
for i = 1:3
fDelta_i = @(d,z) 1./Sigma{i}(z) .* (dSigma{1}(d).*IND{u,1}(d)+dSigma{2}(d).*IND{u,2}(d)+dSigma{3}(d).*IND{u,3}(d)); % Only because we are going to integrate over d from 0 to z
Ups_i = @(z_) arrayfun (@(z) integral(@(d) fDelta_i(d,z) .* (1-P_Ma{M}(d)), zero,z,'ArrayValued',1), z_) ;
dUps_i = @(z_) arrayfun(@(z) fDelta_i(z,z).*(1-P_Ma{M}(z)) - dSigma{i}(z)./(Sigma{i}(z).^2) .* ...
               integral(@(d) dSigma{1}(d).*IND{u,1}(d)+dSigma{2}(d).*IND{u,2}(d)+dSigma{3}(d).*IND{u,3}(d), 0,z,'ArrayValued',1), z_)  ;
ps = @(z) Sigma{i}(z) / (pi*Rd^2) ;
dps = @(z) dSigma{i}(z) / (pi*Rd^2) ;
P_k = @(z,k) nchoosek(NA,k) .* ps(z).^k .* (1-ps(z)).^(NA-k) ;
dP_k = @(z,k) nchoosek(NA,k) .* ps(z).^(k-1) .* (1-ps(z)).^(NA-k-1) .* dps(z) ;
Fbar_Ma_i{u,i} = @(z_) arrayfun(@(z) 0, z_) ;
% f_Ma_i{u,i} = Fbar_Ma_i{u,i} ;
    for k = 0:NA  
         Fbar_Ma_i{u,i} = @(z_) arrayfun(@(z) Fbar_Ma_i{u,i}(z) + nchoosek(NA,k) .* (ps(z)).^k .* (1-ps(z)).^(NA-k) .* (Ups_i(z)).^k, z_) ;   
         F_Ma_i{u,i} = @(z_) arrayfun(@(z) 1-Fbar_Ma_i{u,i}(z), z_) ;
%          f_Ma_i{u,i} = @(z_) arrayfun(@(z) f_Ma_i{u,i}(z) + dP_k(z,k) .* Ups_i(z) + P_k(z,k) .* dUps_i(z) , z_) ; 
    end
    syms z
    f_Ma_i{u,i} = eval( ['@(z)' char( diff( Ups_i(z) ) )] ) 
end
   F_Ma{u,M} = @(z_) arrayfun(@(z) 1- ( min(realmax,Fbar_Ma_i{u,1}(z)) .* IND{u,1}(z) + min(realmax,Fbar_Ma_i{u,2}(z)) .* IND{u,2}(z) ...
                        + min(realmax,Fbar_Ma_i{u,3}(z)) .* IND{u,3}(z) ), z_) ; 
 f_Ma{u,M} = @(z_) arrayfun(@(z) - ( min(realmax,f_Ma_i{u,1}(z)) .* IND{u,1}(z) + min(realmax,f_Ma_i{u,2}(z)) .* IND{u,2}(z) ...
                         + min(realmax,f_Ma_i{u,3}(z)) .* IND{u,3}(z) ), z_) ;
end
end

The problem is that either when I computed manually and wrote it on MATLAB, when I used the symbolic variable, or when I tried to compute the derivative for each one of the pieces, I got errors like 'Index in position 1 exceeds array bounds (must not exceed 2)' or 'A and B must be floating-point scalar' for the function integral (although I'm using arrayfun).

Please let me know any suggested approach to compute this derivative or fix the errors I got.

Thank you in advance!