Implement Variable Step Size for Nonlinear Function MATLAB

I have the following function that returns a nonlinear signal as a function of time:

function y = NARMA_optimized(u,mu,delta,phi,psi,m,N,k)
y = zeros(N,1);
% y(0)=0, and also y(-1)=0 etc.
t = 1;
y(t) = mu.* u(t-m-k).* u(t-k) + delta;
for t = 2:N
    v = 1;
    for i = 1:m
        if t-i > 0  % else we add 0
            v = v + y(t-i);
        end
    end
    v = phi.*y(t-1) + v.*psi.*y(t-1) + mu.*u(t-m-k).*u(t-k) + delta;
    y(t) = v;
end  
end

To call the function, I do the following:

% NARMA Model Comparison
alpha = 0.2;
beta = 0.15;
f_0_i = 1; 
delta = 0.1;
Mu = 1.5;
m = 10;
phi = 0.3;
psi = 0.05;
k_shift = 0;
Sigma = 1;
M = 100; % Upper simulation time limit [seconds]
f = @(t) Sigma.*f_0_i.*sin(alpha.*t).*cos(beta.*t);
dt = 10^-2; % timescale step size [seconds]
t_span = [0:dt:M];
T = NARMA_optimized(f,Mu,delta,phi,psi,m,length(t_span),k_shift);
figure; 
plot(T);

I would like to change the timescale of the current function to match my discretized representation of time denoted as t_span.

I know that this will require two key changes in the function:

1. start at t=0, not t=1.
2. increase the update of the function output y by a step size equivalent to dt, not 1 second. I'm not sure if the current for loop structure can be modified to account for these necessary changes or if I need to start from scratch. I know MATLAB starts at index 1, not index 0. A workaround for "starting" at index 0 in loops is the following:

dt = 10^-2; % [seconds]
M = 100; % [seconds]
t_span = [0:dt:M];
for t=1:numel(t_span)
y(t) = [insert value here]
end

I'm not sure if this is they best way to account for dt and starting at index 0 in the new function.

You can pass t_span into the function directly. The real time will be used for u, or f, if you will. Then define an index for the for loop; this index will be used for y.

function y = NARMA_optimized(u,mu,delta,phi,psi,m,t_span,k)
y = zeros(length(t_span),1);
% y(0)=0, and also y(-1)=0 etc.
%t = 1;
t = t_span(1);
index = 1;
y(index) = mu.* u(t-m-k).* u(t-k) + delta;
for t = t_span(2:end)
    index = index + 1;
    v = 1;
    for i = 1:m
        if index-i > 0  % else we add 0
            v = v + y(index-i);
        end
    end
    v = phi.*y(index-1) + v.*psi.*y(index-1) + mu.*u(t-m-k).*u(t-k) + delta;
    y(index) = v;
end  
end

Change the following in the main code:

M = 1; % Upper simulation time limit [seconds]
T = NARMA_optimized(f,Mu,delta,phi,psi,m,t_span,k_shift);
plot(t_span,T);    

Here is the output: