Index exceeds the number of array elements (0)

clear 
clc
i=100000;
 N=zeros(i:1);
 s=zeros(i:1);
 alfa=zeros(i:1);
 Vs=zeros(i:1);
 k2=zeros(i:1);
 Va=zeros(i:1);
 Vr=zeros(i:1);
 t=zeros(i:1);
 
%parametros


kB=1.36e-25; %constante dealfa (i:1); boltzmann atm L/K
T=291;       % temperatura absoluta ambiente de quito K
u=1.18e-05;  % viscosidad del suero de leche atm*s
E= 99.23; %constante dielectrica
a=500; %radio de micelio de caseina A
Yo=0.012; %potencial de la superficie Volts
kr=0.1; %es el reciproco de espesor de la doble capa A-1
H=1e-22; %constante de Hamaker para micelas J
ko2=4*kB*T/(3*u); % constante cinetica aparente

deltat=0.01;
N(1:1)=10000000;
s(1:1)=2;
t(1:1)=0;
Vs(1:1)=Vr(1:1)+Va(1:1);
Vr(1:1)=((E*a*Yo^2)/2)*log(1+ exp(-kr*a*(s(1:1)-2)));
Va(1:1)=(-H/6)*((2/(S(1:1)^2-4))+(2/S(1:1)^2) + log((S(1:1)^2-4)/s(1:1)^2));
for k=1:i
    funcionN=-k2(k)*(N(k)^2);  %ec1
    N(k+1)= N(k)+ deltat*funcionN; 
 k2(k+1)= ko2/(2*alfa(k+1));      %ec2
   
   funcionalfa=exp((Vs(k)/(kB*T))/s(k)^2); %ec3
   alfa(k+1)=alfa(k)+deltat*funcionalfa;
   Vs(k+1)=Vr(k+1)+Va(k+1);    %ec4
   Vr(k+1)=((E*a*Yo^2)/2)*log(1+ exp(-kr*a*(s(k+1)-2))); 
   Va(k+1)=(-H/6)*((2/(S(k+1)^2-4))+(2/S(k+1)^2) + log((S(k+1)^2-4)/s(k+1)^2));
    
   t(k+1)=t(k)+deltat;    
end
%graficar los resultados
figure(1) % prepara para crear algo y  y x
plot(t,N, 'LineWidth',1.5) % line width da el tamaño de linea
xlabel('tiempo en h') % la etiqueta del eje x
ylabel ('particulas independientes')
grid on % permite cuadricular el grafico
hold on %asi entiende que viene otra figura

Answer 1

[the cyclist]

Everywhere that you wrote

zeros(i:1)

I think you meant

zeros(i,1)