why does matlab report "??? Index exceeds matrix dimensions"?

I have worked on my problem for a long time,but I could not find why matlab report "??? Index exceeds matrix dimensions". Can anybody help me find where index exceeds matrix dimensions,and how to find it. I used the try-catch , but still doesnt work. here is the error message ,"??? Index exceeds matrix dimensions.

Error in ==> quad at 79
if ~isfinite(y(7))

Error in ==> G_yi at 19
f(i)=quad(@sumin,Friction1,Friction2,[],[],i);

Error in ==> quad at 71
y = f(x, varargin{:});

Error in ==> double_int at 11
SS=quad(@G_yi,Protrusion1,Protrusion2,[],[],Friction1,Friction2);  "

I appreciate if anyone can help fix this problem. do I use the try-catch in right way? thank you very much!

here is the code. it includes one main code and 3 function codes.

        Di=4e-3; 
  PR=double_int(0,Di,0,Di) 

function SS=double_int(innlow,innhi,outlow,outhi) 
Protrusion1=outlow;Protrusion2=outhi;Friction1=innlow;Friction2=innhi; 
SS=quad(@G_yi,Protrusion1,Protrusion2,[],[],Friction1,Friction2);

function f=G_yi(Protrusion,Friction1,Friction2) 
Protrusion=Protrusion(:);n=length(Protrusion); 
% if isnumeric(Friction1)==1;FrictionF1=Friction1*ones(size(Protrusion)); 
% else FrictionF1=feval(Friction1,Protrusion); 
% end 
% if isnumeric(Friction2)==1;FrictionF2=Friction2*ones(size(Protrusion)); 
% else FrictionF2=feval(Friction2,Protrusion); 
% end 
save G_yi.mat; 
for i=1:n 
f(i)=quad(@sumin,Friction1,Friction2,[],[],i); 
end 
f=f(:); 

function fun=sumin(Friction,i) 
% global Protrusion; 
% Protrusion(i)=evalin(Protrusion(i)); 
load G_yi.mat; 
MeanShearStress=5 ; %unit Pa 
Density=1e3; 
Ustar=sqrt(MeanShearStress/Density); 
N=15; 
D=zeros(N,1); 
p=zeros(N,1); 
for k=1:N-1; 
D(1)=0.5e-3; 
p(1)=1/N; 
D(k+1)=D(k)+0.5e-3; 
p(k+1)=p(k); 
end; 
%p1=1/3, p2=1/3,p3=1/3 
%D1=40e-6,D2=60e-6, D3=80e-6 
KinematicViscosity=1.004e-6; 
D50=4e-3; 
Roughness=2*D50; 
RoughnessRenolds=Ustar*Roughness/KinematicViscosity; 
if RoughnessRenolds>100 
Intensity=Ustar*2.14; 
Skewness=0.43; 
Flatness=2.88; 
else 
Intensity=Ustar*(-0.187*log(RoughnessRenolds)+2.93); 
Skewness=0.102*log(RoughnessRenolds); 
Flatness=0.136*log(RoughnessRenolds)+2.30; 
end 
CoefficientC=-0.993*log(RoughnessRenolds)+12.36; 
SpecificWeightofSand=1.8836e4; 
SpecificWeightofWater=9.789e3 ; 
Di=4e-3; 
% Protrusion=1e-3; 
Thickness=1.5*D50; 
Y1=0.25*Thickness; 
Y2(i)=0.25*Thickness+Protrusion(i); 
Sum=zeros(N,1); 
for m=1:N-1 
% syms y; 
Kapa=0.4; 
ff1=@(y) (Ustar*CoefficientC.*y/Thickness).*sqrt((0.5*Di)^2-(y-Protrusion(i)-Y1+0.5*Di).^2); 
ff2=@(y) sqrt((0.5*Di)^2-(y-Protrusion(i)-Y1+0.5*Di).^2); 
if Y2(i)<=Thickness 
try 
        MeanBedVelocity(i)=quad(ff1,Y1,Y2(i))./(quad(ff2,Y1,Y2(i))); 

catch
    MeanBedVelocity(i)=quad(ff1,Y1,Y2(i))./(quad(ff2,Y1,Y2(i))+1e-6);
end 
lasterr;else 
MeanBedVelocity(i)=(quad(ff1,Y1,Thickness)+quad(ff1,Thickness,Y2(i)))./(quad(ff2,Y1,Y2(i))); 
end if MeanBedVelocity(i)<=Ustar*CoefficientC 
Yb(i)=(MeanBedVelocity(i)*Thickness)/(Ustar*CoefficientC); 
else 
Yb(i)=Thickness*exp(Kapa*(MeanBedVelocity(i)/Ustar-CoefficientC)); 
end; ParticleRenolds(i)=MeanBedVelocity(i).*Protrusion(i)/KinematicViscosity; 
if ParticleRenolds(i)<=1754 
Cd(i)=(24./ParticleRenolds(i)).*(1+0.15*ParticleRenolds(i).^0.687); 
else 
Cd(i)=0.36; 
end; 
(Protrusion,Friction,Dk,MeanBedVelocity,Yb,Cd, Cl,SpecificWeightofSand,SpecificWeightofWater,MeanShearStress,Ustar,RoughnessRenolds,N,Y1,Di) 
h1(i)=Yb(i)-Y1-Protrusion(i)+0.5*Di; 
h2=Di*(Friction+0.5*D(m)-0.5*Di)/(Di+D(m)); 
Ld=h1(i)+h2; 
Ll=sqrt((0.5*Di)^2-h2.^2); 
Lw=Ll; 
Pei=zeros(N,1); 
Phi=zeros(N,1); 
for j=2:N; 
Pei(1)=p(1)*Di/(Di+D(1)); 
Pei(j)=p(j)*Di/(Di+D(j))+ Pei(j-1); 
Phi(j)=1-Pei(j); 
end; 
HidingFactor=(Pei(N)/Phi(N))^0.6; 
EffectiveShearStress=HidingFactor*MeanShearStress; 
% SpecificWeightofSand=1.8836e4; 
% SpecificWeightofWater=9.789e3 ; %at 20 degree centigrade 
DimensionlessEffectiveShearStress= EffectiveShearStress/((SpecificWeightofSand-SpecificWeightofWater)*Di); 
ff3=@(y) sqrt((0.5*Di)^2-(y-Protrusion(i)-Y1+0.5*Di).^2); 
         try
          A(i)=quad(ff3,Y1,Y2(i));
         catch
            A(i)=A(i)+1e-6  ;
         endBr=Ustar*sqrt((2*Lw*pi*Di^2)./((Cd(i).*Ld+Cl(i).*Ll)*6.*A(i)*DimensionlessEffectiveShearStress)); % Rolling Threshold 
Bl=Ustar*sqrt((2*pi*Di^2)/(Cl(i)*6.*A(i)*DimensionlessEffectiveShearStress)); %Lifting Threshold %syms Ub ; % Instantaneous velocity 
% U=((Ub-MeanBedVelocity)/Intensity);% PDF=@(Ub) (exp(-((Ub-MeanBedVelocity)/Intensity).^2/2)/(Intensity*sqrt(2*pi))).*(1+(Skewness/factorial(3))*(((Ub-MeanBedVelocity)/Intensity).^3-3*((Ub-MeanBedVelocity)/Intensity))+(Flatness-3)*(((Ub-MeanBedVelocity)/Intensity).^4-6*((Ub-MeanBedVelocity)/Intensity).^2+3)/factorial(4))/Intensity; %PDF of Velocity Fluctuation 
% Pr(i)=quad(PDF,-Bl,-Br)+quad(PDF,Br,Bl); 
% Pr(i)=sum(arrayfun(@(BL,BR) quad(PDF, -BL, -BR) + quad(PDF, BR, BL), Bl, Br)); 
syms Ub ; % Instantaneous velocity 
  PD(i)=int((exp(-((Ub-MeanBedVelocity(i))/Intensity).^2/2)/(Intensity*sqrt(2*pi))).*(1+(Skewness/factorial(3))*(((Ub-MeanBedVelocity(i))/Intensity).^3-3*((Ub-MeanBedVelocity(i))/Intensity))+(Flatness-3)*(((Ub-MeanBedVelocity(i))/Intensity).^4-6*((Ub-MeanBedVelocity(i))/Intensity).^2+3)/factorial(4))/Intensity,-Bl,-Br)...
               +int((exp(-((Ub-MeanBedVelocity(i))/Intensity).^2/2)/(Intensity*sqrt(2*pi))).*(1+(Skewness/factorial(3))*(((Ub-MeanBedVelocity(i))/Intensity).^3-3*((Ub-MeanBedVelocity(i))/Intensity))+(Flatness-3)*(((Ub-MeanBedVelocity(i))/Intensity).^4-6*((Ub-MeanBedVelocity(i))/Intensity).^2+3)/factorial(4))/Intensity,Br,Bl);Pr(i)= double(PD(i)); 
% Pl=1-quad(PDF,-Bl,Bl); %end Sum(1)=p(1)*Pr(i); 
Sum(m+1)=p(m)*Pr(i)+Sum(m); 
end; 
fun=Sum(m+1);

[Grzegorz Knor]

The code is completely unreadable and impossible to copy.

http://www.mathworks.com/matlabcentral/answers/help/markup

[Lin LI]

thank you for your suggestions, I've changed it. I didnt know how to do it before. thank you.

I appreciate very much if you can help fix this problem.

Answer 1

[Walter Roberson]

Read the documentation for quad():

The function y = fun(x) should accept a vector argument x and return a vector result y, the integrand evaluated at each element of x.

Does that match your function sumin() ?