This is the code thus far: (All values defined, just not shown)

t = 0:dt:tmax; % time
NK = zeros(size(t)); % number of unphosphorylated kinase
NKP = zeros(size(t)); % number of phosphorylated kinase
NP = zeros(size(t)); % number of phosphatase

% Now put intial values into first element of the vectors
NK(1) = NK_init;
NKP(1) = NK_total-NK(1);
NP(1) = NP_init;

% Now loop through time
for i = 2:length(t)

rphosph = NK(i-1)*(kbase + kauto*NKP(i-1)); % total phosphorylation rate

rdephos = NKP(i-1)*k_dephos*NP(i-1)/(KM + NKP(i-1)); % total dephosphorylation rate

NKP(i) = NKP(i-1) + dt*(rphosph-rdephos); % update number of phosphorylated kinase
NK(i) = NK_total - NKP(i); % update number of dephosphorylated kinase
NP(i) = NP(i-1); % number of phosphatase unchanged

end

How would I introduce an outer loop in the code, so that the initial number of unphosphorylated
kinases steps from the maximum number (50) down to zero?