Cyclic Voltammetry Simmulation in MATLAB, I am running into issues with my data points returning as NaN values, i am a beginner, any help wanted

Question

I am having a broad degree of issues resulting in NaN values for a large amount of my output file.

I have tried changing different variables one by one without sucess and have yet to find the root of the issues. When the output file reads the x values representing potential are present but the current are all NaN, i have broadened the search for mistaken variables by having the variable values print during loops, but it has just confused me more, as nearly all of the values are NaN. I appreciate any help possible thank you.

My code is shown below:

% Cyclic voltammetry simulation


% Run the input script
run('inputs.m');


% Concentration of all species (mol/dm**3)
c = [5.200000E-03, 0, 0];

% Set Cell Resistance
cellr = 10000;

% Diffusion Coefficients of all Species (cm**2/s)
d = [1.000000E-05, 1.000000E-05, 1.000000E-05];

% Scan Rate (V/s)
sr = 0.1;

% Temperature (K)
tk = 295;

% Heterogeneous Rate Constant (cm/s)
rhet1 = 0.5;

% Homogeneous Rate Constant (/M/s)
rhom = 0;

% Area of Electrode (cm**2)
ae = 0.07065;

% Initial and Final Potentials (V)
pi = 0.685;
pf = -0.915;

% Transfer Coefficients
al1 = 0.5;

% Background Slope(SL) (mA/V) and Y-intercept (YI) (mA)
sl = 0;
yi = 0;

% Number of Increments before Reversal
mr = 5000;

% Simulator Diffusion Coefficients
dx = 0.45;


% Enter file name
a = input('Enter file name: ', 's');
fid = fopen(a, 'w');

% Initialize bg
bg = 0;

% The dimensioning is staying at 4 ñ previous program
dd = zeros(1,3);
zz = zeros(1,3);
d = zeros(1,3);
c = zeros(1,3);
u2 = zeros(3,5500);
u3 = zeros(3,5500);

% Simulator Diffusion Coefficients
dx = 0.45;

% Simulation Constants
% ********Subroutine to identify DMAX *************
dmax = max(d);

% ********* Time (TT) and Distance (X) Units ******
tt = abs(pi-pf)/(mr*sr);
xx = sqrt((tt*dmax)/dx);

% Total Number of Increments 
ni = mr*2;

% Number of Species
ns = 3;

const = 96486.7/(8.314*tk);
sp = (pf-pi)/mr;
nm = 6*sqrt(dx*ni);

% ********** SET INITIAL CONDITIONS ***************
% Potential
pv = pi;

% Iteration Number
nt = 1;

% Maximum Current 
zmax = 0;

% Diffusion Coefficients
dd(1) = dx*d(1)/dmax;
dd(2) = dx*d(2)/dmax;
dd(3) = dx*d(3)/dmax;

% Heterogenous Rates
rh1 = (rhet1*tt/xx);

% Homogenous Rates
hom = rhom*tt;

% Initial Concentrations
for j = 1:ns
    for i = 1:nm+10
        u3(j,i) = c(j);
        u2(j,i) = c(j);
    end
end

% ************************************
%         START OF TIME LOOP
% ************************************
while nt ~= (ni+1)
    % Calculate Fluxes at Electrode Surface
    if (nt == (mr+2))
        sp = -sp;
    end
    if (nt > 1)
        pv = pv+sp;
    end
    rf1 = rh1*exp(-al1*pv*const);
    rb1 = rh1*exp((1-al1)*pv*const);
    aa = rf1*u2(1,1);
    bb = rb1*u2(2,1);
    qq = rb1/(2*dd(2));
    rr = rf1/(2*dd(1));
    zz(1) = ((aa-bb)/(1+rr+qq));
    zz(2) = -zz(1);
    zf = zz(1);
    zr = (zf*96486.7*ae*sqrt(sr*dmax)/(sqrt(dx*abs(sp))));
    % CURR is current in mA
    curr = zr+bg;
    curra = curr/1000;
    % shift voltage based on IR drop
    shift = curra*cellr;
    newpv = pv-shift;
    % CURRA is current in A
    fprintf(fid, '%f %f\n', pv, curra);

       % Print the values of the variables
    fprintf('rf1: %f, rb1: %f, u2(1,1): %f, u2(2,1): %f, dd(1): %f, dd(2): %f\n', rf1, rb1, u2(1,1), u2(2,1), dd(1), dd(2));
   
    % Calculate Background
    if nt == mr
        yi = -yi;
    end
    bg = (sl*(-pv)+yi);

    % DIFFUSION BLOCK
    nm = 6*sqrt(dx*nt);
    for j = 1:ns
        u3(j,1) = u2(j,1)+dd(j)*(u2(j,2)-u2(j,1))-zz(j);
        for i = 2:nm
            u3(j,i) = u2(j,i)+dd(j)*(u2(j,i+1)-2*u2(j,i)+u2(j,i-1));
        end
    end
    for i = 1:nm
        u2(1,i) = u3(1,i);
        u2(2,i) = u3(2,i);
        u2(3,i) = u3(3,i);
    end

    % PRINT NT,NI
    nt = nt+1;
end
fclose(fid);

figure; % Create a new figure window
plot(pv, curr); % Plot current as a function of potential
title('Current vs Potential'); % Add a title to the plot
xlabel('Potential (V)'); % Label the x-axis
ylabel('Current (mA)'); % Label the y-axis
grid on; % Add a grid to the plot