循坏索引超出数组范围

%输入 b:坡度 h:坡高 gama：重度 phi0:内摩擦角 c:粘聚力
%输出 Fsmin x0 y0 R0
Fs=2.0;
b=0.5;
h=25;
gama=20;
phi0=19.6;
c=3;
phi0=phi0*pi/180;
alpha=atan(b);
L=h/sin(alpha);%坡长
m=L*cos(alpha);
beta=[500,500];
x=[];
y=[];
xd=[];
yd=[];
qs=[];
T=[];
P=[];
l=[];
P(1)=0;
y2=[];
h1=[];
W1=[];
for x=-2*h:0.1:h %圆心横坐标取值范围
    for y=h:0.1:3*h %圆心纵坐标取值范围
        R=sqrt(x^2+y^2); %新圆半径
        x1=sqrt(R^2-(h-y)^2)+abs(x); %x最大取值
        d1=x1/500; %每块土条宽度
       for xd=x1:-d1:0 %圆弧曲线的x坐标取值
       yd=y-sqrt(R^2-(x-xd)^2); %圆弧所在方程
       beta=atan((xd-x)/(y-yd));%土条坡底角度
if xd<=m 
  y2=tan(alpha)*xd; %斜坡方程
  h1=abs(y2-yd); %土条高度
  W1=gama*h1*d1; %土条的重力
else
   h1=abs(h-yd);
   W1=gama*h1*d1; %土条的重力
end
       end
  s=1;
  n=500;
  for i=s:n   
      if i==1
        l(i)=d1/cos(beta(i)); %斜底面长度
        qs(i)=0;
        T(i)=(c*l(i)+W1(i)*cos(beta(i))*tan(phi0)/Fs);%底面剪力
        P(i)=W1(i)*sin(beta(i))-T(i); %条间推力  
        i=i+1;
      else
  l(i)=d1/cos(beta(i)); %斜底面长度
  qs(i)=cos(beta(i-1)-beta(i))-tan(phi0)*sin(beta(i-1)-beta(i))/Fs; %传递系数
  T(i)=(c*l(i)+W1(i)*cos(beta(i))*tan(phi0)/Fs);%底面剪力
  P(i)=W1(i)*sin(beta(i))+qs(i)*P(i-1)-T(i); %条间推力
      end
  i=i+1;
 
  end
    P;
    if P(n)>0.5
      Fs=Fs-0.0001;
    elseif P(n)<-0.5
      Fs=Fs+0.0001;
    else
        disp(Fs);
    end
    end
end

以上是我的代码，报错显示
索引超出数组范围。

出错 nnn6666 (line 53)
l(i)=d1/cos(beta(i)); %斜底面长度