A never used variable's index out of bounds! From the function quadl? weird, dont know why.

I run this code and it gave me error message, and it seems that this problem comes from the function quadl, because the variable y never appear in my code, y is from the function quadl.I set some breakpoint,still couldnt find why.If anyone can help fix this problem , i appreciate it very much. anyone can help find what the problem is. I am really annoyed by this problem.thank you so much!

here is the error:

      "??? Attempted to access y(13); index out of bounds because numel(y)=1.

Error in ==> quadl at 72
if ~isfinite(y(13))

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

Error in ==> quadl at 64
y = feval(f,x,varargin{:}); y = y(:).';

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

Error in ==> main at 4
PR=double_int(0,Di,0,Di) "

here is the code, one main code and 3 functions.

      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=quadl(@G_yi,Protrusion1,Protrusion2,[],[],Friction1,Friction2);  

    function f=G_yi(Protrusion,Friction1,Friction2)
  Protrusion=Protrusion(:);n=length(Protrusion);
  save G_yi.mat;
  f=zeros(n,1);
  for i=1:n
  f(i)=quadl(@sumin,Friction1,Friction2,[],[],i);
  end
  f=f(:);

    function fun=sumin(Friction,i)

    load G_yi.mat;
  MeanShearStress=5 ;   
  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;

    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;

    Thickness=1.5*D50;
      Y1=0.25*Thickness;
      Y2(i)=0.25*Thickness+Protrusion(i);

      Sum=zeros(N,1);
  for m=1:N-1

      Kapa=0.4;
      ff1=@(t) (Ustar*CoefficientC.*t/Thickness).*sqrt((0.5*Di)^2-(t-Protrusion(i)-Y1+0.5*Di).^2);
      ff2=@(t) sqrt((0.5*Di)^2-(t-Protrusion(i)-Y1+0.5*Di).^2);
      MeanBedVelocity=zeros(n,1);
      MM=zeros(n,1);
       MM(i)= quadl(ff2,Y1,Y2(i));
      if Y2(i)<=Thickness

           MeanBedVelocity(i)=quadl(ff1,Y1,Y2(i))./(MM(i)); 

      else
          MeanBedVelocity(i)=(quadl(ff1,Y1,Thickness)+quadl(ff1,Thickness,Y2(i)))./(quadl(ff2,Y1,Y2(i)));
      end

   Yb=zeros(n,1);
      if MeanBedVelocity(i)<=Ustar*CoefficientC
          Yb(i)=(MeanBedVelocity(i)*Thickness)/(Ustar*CoefficientC);
      else
          Yb(i)=Thickness*exp(Kapa*(MeanBedVelocity(i)/Ustar-CoefficientC));
      end;
       ParticleRenolds=zeros(n,1);
      ParticleRenolds(i)=MeanBedVelocity(i).*Protrusion(i)/KinematicViscosity;
       Cd=zeros(n,1);
      if ParticleRenolds(i)<=1754
          Cd(i)=(24./ParticleRenolds(i)).*(1+0.15*ParticleRenolds(i).^0.687);
      else
          Cd(i)=0.36;
      end;
    Cl=zeros(n,1);     

          Cl(i)=Cd(i);

          h1=zeros(n,1);    
  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;

            DimensionlessEffectiveShearStress= EffectiveShearStress/((SpecificWeightofSand-SpecificWeightofWater)*Di); 
            ff3=@(t) sqrt((0.5*Di)^2-(t-Protrusion(i)-Y1+0.5*Di).^2);
          A=zeros(n,1);

           A(i)=quadl(ff3,Y1,Y2(i));

            Br=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  ; 
  Pr=zeros(n,1);
   PD=zeros(n,1);

             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));

               Sum(1)=p(1)*Pr(i);
               Sum(m+1)=p(m)*Pr(i)+Sum(m);
   end;
           fun=Sum(m+1);

Answer 1

[Walter Roberson]

The function passed to quad must take exactly one parameter. Your sumin() takes two parameters.

Does your sumin() return a vector with as many results as there were values in the single input parameter?