Ошибка - Index exceeds the number of array elements. Index must not exceed 1

clc
clear all
%%исходные данные
Vx(1) = -2954; % проекция скорости на ось х
Vy(1) = -759; % проекция скорости на ось у
Vz(1) = 1965; % проекция скорости на ось z
V(1) = sqrt(Vx^2 + Vy^2 + Vz^2); % абсолютная скорость
p(1) = 0.412; %плотность воздуха(h = 10000 м)
x(1) = 1831545; % координата х
y(1) = 5966464; % координата у
z(1) = -1673435; % координата z
g(1) = 9.8; % ускорение свободного падения, м/с^2
Sm = 3.6361; % площадь Миделя,м^2
m = 1000; % масса, кг
Cx = 0.3; % коэффициент лобового сопротивления
R = 6371000; % средний радиус Земли,м
r(1) = sqrt((R + y)^2 + x^2 + z^2); % истинное расстояние точки от центра Земли,м
W = 0.00007292; % угловая скорость Земли, с^-1
k = (Sm/2*m);
fi = 0.87; % широта, рад
fu = 0.7; % долгота, рад
a11 = W^2*(sin(fi)^2+cos(fi)^2*sin(fu)^2);
a12 =-W^2*sin(fi)*cos(fi)*cos(fu);
a13 = W^2*cos(fi)^2*sin(fu)*cos(fu);
a21 = -W^2*sin(fi)*cos(fi)*cos(fu);  
a22 = W^2*cos(fi)^2;
a23 = W^2*sin(fi)*cos(fi)*sin(fu);
a31 = W^2*cos(fi)^2*sin(fu)*cos(fu);
a32 = W^2*sin(fi)*cos(fi)*sin(fu);
a33 = W^2*(sin(fi)^2+cos(fi)^2*cos(fi)^2);
b12 = -2*W*cos(fi)*sin(fu);
b13 = -2*W*sin(fi);
b21 = 2*W*cos(fi)*sin(fu);
b31 = 2*W*sin(fi);
b23 = 2*W*cos(fi)*cos(fu);
b32 = -2*W*cos(fi)*cos(fu);
Vxtilda(1) = Vx(1); % только на первом шаге интегрирования
Vytilda(1) = Vy(1); % только на первом шаге интегрирования
Vztilda(1) = Vz(1); % только на первом шаге интегрирования
dt = 0.1; % шаг интегрирования, с
t(1) = 0; % начальное время, с
i = 1; % счётчик    
FVx(1) = 0;
 
while (r(i) > R)
 
    FVx(1) = 0;
 
    t(i+1) = t(i) + dt;
 
    FVx(i) = -k * Cx * p(i) * V(i) * Vx(i) - (g(i)/r(i)) * x(i) + a11 * x(i) + a12 * (R + y(i)) + a13 * z(i) + b12 * Vy(i) + b13 * Vz(i);
 
    FVxtilda(i) = -k * Cx * p(i) * V(i) * Vxtilda(i) - (g(i)/r(i)) * x(i) + a11 * x(i) + a12 * (R + y(i)) + a13 * z(i) + b12 * Vy(i) + b13 * Vz(i);
 
    Vxtilda(i+1) = Vx(i) + dt * FVx(i); % дальше
 
    Vx(i+1) = Vx(i) + dt/2 * (FVx(i) + FVxtilda(i));
 
    x(i+1) = x(i) + dt * Vx(i);
 
 
    FVy(i) = -k * Cx * p(i) * V(i) * Vy(i) - (g(i)/r(i))*(R + y(i)) + a21 * x(i) + a22 * (R + y(i)) + a23 * z(i) + b21 * Vx(i) + b23 * Vz(i);
 
    FVytilda(i) = -k * Cx * p(i) * V(i) * Vytilda(i) - (g(i)/r(i))*(R + y(i)) + a21 * x(i) + a22 * (R + y(i)) + a23 * z(i) + b21 * Vx(i) + b23 * Vz(i);
 
    Vytilda(i+1) = Vy(i) + dt * FVy(i); % дальше
 
    Vy(i+1) = Vy(i) + dt/2 * (FVy(i) + FVytilda(i));
 
    y(i+1) = y(i) + dt * Vy(i);
 
 
    FVz(i) = -k * Cx * p(i) * V(i) * Vz(i) - (g(i)/r(i)) * z(i) + a31 * x(i) + a32 * (R + y(i)) + a33 * z(i) + b31 * Vx(i) + b32 * Vy(i);
     
    Fvztilda(i) = -k * Cx * p(i) * V(i) * Vztilda(i) - (g(i)/r(i)) * z(i) + a31 * x(i) + a32 * (R + y(i)) + a33 * z(i) + b31 * Vx(i) + b32 * Vy(i);
 
    Vztilda(i+1) = Vz(i) + dt * FVz(i); %дальше
 
    Vz(i+1) = Vz(i) + dt/2 * (FVz(i) + Fvztilda(i));
 
    z(i+1) = z(i) + dt * Vz(i);
    
    r(i+1) = sqrt((R + y(i+1))^2 + x(i+1)^2 + z(i+1)^2);
 
    i = i+1
end