This code calculate m n h values in For loop, and m n h values are probability values, they must be between 0 and 1.I found these values.But I want to find spikes values during 100 second.Spike means values over 20mV in this code.

   clear all;clc;
  % === simulation time === 
  Simulationtime= 2000; % as ms
  deltaT=0.01;
  t=0:deltaT:Simulationtime;

% === specify external current I === 

changeTimes = [0]; % as ms% 

CurrentLevel = [50]; %Change this to see the effect of the different currents on the voltage (Suggested values: 3, 20, 50, 1000)

%Set externally applied current over time

%Here, the first 500 time-outs are in the current 50's, the next 1500 hours

%zero current (zero neuron rest potential)

%and the remaining timeouts are constant current

I(1:500) = CurrentLevel; 

I(501:2000) = 0; 

I(2001:numel(t)) =  CurrentLevel;

%Interpret the line above and subtract the following line for constant current

%and the effect of the voltage time

%===fixed parameters===%

g_K=36;   

g_Na=120; 

g_L=0.3;

E_K =-12;   %-12

E_Na=115;   %115

E_L=10.6;   %10,6

C=1;

%===initial situations===%

V=0; % Starting line voltage

alfa_n = 0.01 * ( (V+55) / (1-exp(-(V+55)/10)) ); 

beta_n = 0.125*exp(-(V+65)/80);

alfa_m = 0.1*( (V+40) / (1-exp(-(V+40)/10)) ); 

beta_m = 4*exp(-(V+65)/20);

alfa_h = 0.07*exp(-(V+65)/20);

beta_h = 1/(exp(-(V+35)/10)+1); 

n(1) = alfa_n/(alfa_n+beta_n);

m(1) = alfa_m/(alfa_m+beta_m); 

h(1) = alfa_h/(alfa_h+beta_h); 

i =1

while i ~= numel(t) 

        %Coefficients, currents and derivative calculations in each step.

        %---calculate coefficients---%

        %The equations in this section are the same as above, only calculated in each step.

        alfa_n(i) = 0.01 * ( (10-V(i)) / (exp((10-V(i))/10)-1) );

        beta_n(i) = 0.125*exp(-V(i)/80);

        alfa_m(i) = 0.1*( (25-V(i)) / (exp((25-V(i))/10)-1) );

        beta_m(i) = 4*exp(-V(i)/18);

        alfa_h(i) = 0.07*exp(-V(i)/20);

        beta_h(i) = 1/(exp((30-V(i))/10)+1);

        %---calculate currents---%

        I_Na = (m(i)^3) * g_Na * h(i) * (V(i)-E_Na); 

        I_K = (n(i)^4) * g_K * (V(i)-E_K); 

        I_L = g_L *(V(i)-E_L); 

        I_iyon = 0 - I_K - I_Na - I_L; 

        %---Euler calculates derivative using first order approximation.---%

        V(i+1) = V(i) + deltaT*I_iyon/C;

        gm=sqrt(((2*alfa_m(i)*beta_m(i))/(60*128*(alfa_m(i)+beta_m(i)))))*randn;

        gn=sqrt(((2*alfa_n(i)*beta_n(i))/(18*128*(alfa_n(i)+beta_n(i)))))*randn;

        gh=sqrt(((2*alfa_h(i)*beta_h(i))/(60*128*(alfa_h(i)+beta_h(i)))))*randn;

        n(i+1) = n(i) + deltaT*(alfa_n(i) *(1-n(i)) - beta_n(i) * n(i))+ gm; 

        m(i+1) = m(i) + deltaT*(alfa_m(i) *(1-m(i)) - beta_m(i) * m(i))+ gn;

        h(i+1) = h(i) + deltaT*(alfa_h(i) *(1-h(i)) - beta_h(i) * h(i))+ gh; 

%       disp(n(i+1));

%       disp(m(i+1));

%       disp(h(i+1));

         if n(i+1) < 1 && n(i+1) > 0 %%n değerini istenilen aralıkta tutabilirsin

          if m(i+1) < 1 && m(i+1) > 0  %%m değerini istenilen aralıkta tutabilirsin

            if h(i+1) < 1 && h(i+1) > 0 %%h değerini istenilen aralıkta tutabilirsin

             i = i+1;

            end

          end

         end

end

%%Max ve Min values

max(m)

   max(h)

   max(n)

   min(m)

   min(n)

   min(m)
E_Na_mu = 115;

E_Na_std = 5;

E_Na = normrnd(E_Na_mu, E_Na_std);

V = V-70; %Set resting potential to -70mv

figure

%===plot Voltage===%

plot(t,V,'LineWidth',3)

%hold on

legend({'Voltaj'})

ylabel('Voltaj (mv)')

xlabel('Zaman (ms)')

title('Zamanla Gerilim')

%===plot Conductance===%

%figure

%p1 = plot(t,g_K*n.^4,'LineWidth',2);

%hold on

%p2 = plot(t,g_Na*(m.^3).*h,'r','LineWidth',2);

%legend([p1, p2], 'Potasyum İçin İletkenlik', 'Sodyum İçin İletkenlik')

%ylabel('İletkenlik')

%xlabel('Zaman (ms)')

%title('Potasyum ve Sodyum İyonları İçin İletkenlik')

[ayten kurnaz]

And then I want to save the spikes .dat file.Can you help me ?