Index exceeds the number of array elements error

In my code I am trying to graph a variable magnetic field and I use the rk2 method to solve, but it throws the following error:

Index exceeds the number of array elements (100001).

Error in programa_reto (line 52)
while t(i) < tf

My code is:

%Programa del reto
%Tenemos que F=qvxB, ma=qvxB, a=qvxB/m, d^2r/dt^2=q/m*vxB, donde r es el
%vector de posición.
%Establecemos carga, masa, tiempo inicial, tiempo final y h:
clear, clc, close all

graph2D = false; %Cambiar "False" para graficar en 3D
q = 0.5; m = 0.1; t0 = 0; tf = 10; h = 0.0001; %Cambiar valores

%Inicializando vectores en cero.
N = round((tf - t0) / h); %Define el tamaño de los vectores 

%Vectores de posición
rx = zeros(1, N + 1); 
ry = zeros(1, N + 1); 
rz = zeros(1, N + 1);

%Vectores de posición
vx = zeros(1, N + 1); 
vy = zeros(1, N + 1); 
vz = zeros(1, N + 1);

%Campo magnético
bx = zeros(1, N + 1); 
by = zeros(1, N + 1); 
bz = zeros(1, N + 1);

%Tiempo
t = zeros(1, N + 1);

%Ecuaciones diferenciales
%r"=q/m*vxB=[q/m*(vy*Bz-vxBz)i+(vz*Bx-vx*Bz)j+(vx*By-vy*Bx)k]
%rx"=q/m*(vy*Bz-vxBz)   ry"=q/m*(vz*Bx-vx*Bz)   rz"=q/m*(vx*By-vy*Bx)
%Cambio de variable r´= v
%rx´=vx       ry´=vy         rz´=vz
%vx´= q/m*(vy*Bz-vxBz)   vy´= q/m*(vz*Bx-vx*Bz)   vz´= q/m*(vx*By-vy*Bx)
%Declarando la función para cada componente de v´
f=@(va, Ba, vb, Bb) q * (va * Ba - vb * Bb)/m;
%Declararndo los componentes del campo magnético B=(Bx, By, Bz)

%Nota cambiar según nuestras funciones
Bx = @(x, y, z) sin(-y);
By = @(x, y, z) cos(x^2);
Bz = @(x, y, z) exp(sin(x-y + z));

%Declarar los valores iniciales rx, ry, rz, vx, vy, vz cuando t=t0:
rx(1) = 0; ry(1) = 0; rz(1) = 0; vx(1) = 0.1; vy(1) = 0.2; vz(1) = 0.1; 
t(1) = t0;

%Método RK2
i = 1;
while t(i) < tf
    %Constante k1, m1 para x, y, z
    k1x = vx(i);
    k1y = vy(i);
    k1z = vz(i);
    
    %Constante K1 para vx, vy, vz
    bx(i) = Bx(rx(i), ry(i), rz(i));
    by(i) = By(rx(i), ry(i), rz(i));
    bz(i) = Bz(rx(i), ry(i), rz(i));
    
    k1vx = f(vy(i), bz(i), vz(i), by(i));
    k1vy = f(vz(i), bx(i), vx(i), bz(i));
    k1vz = f(vx(i), by(i), vy(i), bx(i));
    
    %Constante K2
    %m2 para x, y, z
    k2x = vx(i) + h * k1vx;
    k2y = vy(i) + h * k1vy;
    k2z = vz(i) + h * k1vz;
    
    %K2 para vx, vy, vz
    bx(i) = Bx(rx(i) + h*k1x, ry(i) + h*k1y, rz(i) + h*k1z);
    by(i) = By(rx(i) + h*k1x, ry(i) + h*k1y, rz(i) + h*k1z);
    bz(i) = Bz(rx(i) + h*k1x, ry(i) + h*k1y, rz(i) + h*k1z);
    
    k2vx = f(vy(i) + h * k1vy, bz(i), vz(i) + h * k1vz, by(i));
    k2vy = f(vz(i) + h * k1vz, bx(i), vx(i) + h * k1vx, bz(i));
    k2vz = f(vx(i) + h * k1vx, by(i), vy(i) + h * k1vy, bx(i));
    
    rx(i + 1) = rx(i) + h * (k1x + k2x)/2;
    ry(i + 1) = ry(i) + h * (k1y + k2y)/2;
    rz(i + 1) = rz(i) + h * (k1z + k2z)/2;
    vx(i + 1) = vx(i) + h * (k1vx + k2vx)/2;
    vy(i + 1) = vy(i) + h * (k1vy + k2vy)/2;
    vz(i + 1) = vz(i) + h * (k1vz + k2vz)/2;
    
    t(i + 1) = t(i) + h;
    i = i + i; 
end
%Gráfica del campo mgnético (se tienen que recortar los vectores de cada
%"n" elementos para que la gráfica no se sature de vectore)
%axis tight

partes = 5; %Cada arreglo se dividirá en partes, entre más grande sea el 
%valos más vectores aparecerán.
N
n = round(N / partes)
%Se obtienen arreglos muestreados cada "n" elementos. "Rec" es de
%"recortado".
bxRec = bx(1 : n : end); byRec = by(1 : n : end); bzRec = bz(1 : n : end);
rxRec = rx(1 : n : end); ryRec = ry(1 : n : end); rzRec = rz(1 : n : end);

scaleB = 0.2; %Tamaño de los vectores del campo magnético
%Tamaño de los vectores del campo magnético
if graph2D == false
    quiver3(rxRec, ryRec, rzRec, bxRec, byRec, bzRec,scaleB, "r")
else
    quiver(rxRec, ryRec, bxRec, byRec,scaleB, "r")
end
hold on

%Gráfica de la trayectoria
grosor = 1;
if graph2D == false
    plot3(rx, ry, rz, "k", "LineWidth", grosor)
else
    plot(rx, ry, "k", "LineWidth", grosor)
end

legend("B", "r")
xlabel("x")
ylabel("y")
zlabel("z")

Answer 1

[Ameer Hamza]

Change the

i = i + i;

to

i = i + 1;