Infinite Loop, Values not Updating

In the while loop, I am comparing the accuracy of Euler's and Midpoint methods (the methods themselves should be correct). If they differ by a certain amount, I reduce the step size by half and redo the calculation for that time point. Else, I update the indices and move on. Plot(t, x(:,4), t, y(:,4)) should gnerate a graph that portrays an action potential.

Somehow after dif hits the specified threshold, x, y, and dif stop to update. The step size keeps getting reduced and I get stuck in an infinite loop. Does anyone have any insights on what might have happened there? Any thought is appreciated!

dt = 0.08;
Vm = 0;

t(1)=0;

i_inj = 17;

%initial values for n, h, m, V
x = [0.317676914060697	0.596120753508460	0.0529324852572496	0]; %set first value to steady state to 0mV
y = [0.317676914060697	0.596120753508460	0.0529324852572496	0];
i = 1;

k = 1;

while t(k)<11
    
    %Euler's method
    x(i+1,:) = x(i,:)+dt.*derivs(t(k),x(i,:),i_inj);

    %midpoint method
    k1 = dt.*derivs(t(k),y(i,:),i_inj);
    k2 = dt.*derivs(t(k)+0.5*dt,y(i,:)+0.5.*k1,i_inj);
    y(i+1,:) = k2+y(i,:);
    
    dif = abs(x(i+1,4)-y(i+1,4));
    
    if (dif>=0.1)
        dt = dt*0.5;
        
    else
        t(k+1) = dt+t(k);
        i = i+1;
        k = k+1;
    end 
    
end


function y = derivs(t,x,I_inj)
g_Kbar = 36;
E_K = -12;
g_Nabar = 120;
E_Na = 115;
Gm = 0.3;
V_rest = 10.613;
V=x(4);


if(V == 10)
    alpha_n = 1/10;
else
   alpha_n = (10-V)./(100*(exp(1-V/10)-1)); 
end
beta_n = 0.125*exp(-V/80);
n=x(1);
dndt = alpha_n*(1-n)-beta_n*n;
y(1) = dndt; 


alpha_h = 0.07*exp(-V/20);
beta_h = 1/(exp(3-V/10)+1);
h=x(2);
dhdt = alpha_h*(1-h)-beta_h*h;
y(2) = dhdt; 


if (V==25)
    alpha_m = 1;
else
    alpha_m = (25-V)/(10*(exp(2.5-V/10)-1));
end
beta_m = 4*exp(-V/18);
m=x(3);
dmdt = alpha_m*(1-m)-beta_m*m;
y(3) = dmdt; 


if t>=1 && t<=1.5
    dVdt = I_inj+g_Nabar.*m.^3.*h.*(E_Na-V)+g_Kbar.*n.^4.*(E_K-V)+Gm*(V_rest-V);
else
    dVdt = g_Nabar.*m.^3.*h.*(E_Na-V)+g_Kbar.*n.^4.*(E_K-V)+Gm*(V_rest-V);
end

y(4) = dVdt;
end

[KSSV]

What is this function derivs?

[Irene Zhou]

It's a derivative function that calculates and spits out four variables in a row. This function should be correct as I have tested it in other places.

function y = derivs(t,x,I_inj)
g_Kbar = 36;
E_K = -12;
g_Nabar = 120;
E_Na = 115;
Gm = 0.3;
V_rest = 10.613;
V=x(4);


if(V == 10)
    alpha_n = 1/10;
else
   alpha_n = (10-V)./(100*(exp(1-V/10)-1)); 
end
beta_n = 0.125*exp(-V/80);
n=x(1);
dndt = alpha_n*(1-n)-beta_n*n;
y(1) = dndt; 


alpha_h = 0.07*exp(-V/20);
beta_h = 1/(exp(3-V/10)+1);
h=x(2);
dhdt = alpha_h*(1-h)-beta_h*h;
y(2) = dhdt; 


if (V==25)
    alpha_m = 1;
else
    alpha_m = (25-V)/(10*(exp(2.5-V/10)-1));
end
beta_m = 4*exp(-V/18);
m=x(3);
dmdt = alpha_m*(1-m)-beta_m*m;
y(3) = dmdt; 


if t>=1 && t<=1.5
    dVdt = I_inj+g_Nabar.*m.^3.*h.*(E_Na-V)+g_Kbar.*n.^4.*(E_K-V)+Gm*(V_rest-V);
else
    dVdt = g_Nabar.*m.^3.*h.*(E_Na-V)+g_Kbar.*n.^4.*(E_K-V)+Gm*(V_rest-V);
end

y(4) = dVdt;
end

Answer 1

[Andy]

I don't have real values for n0 etc so I put some random numbers in and the resultant t(k) was generally small, much less than the threshold of 11 that you have set for the while loop. I was just wondering if this is the correct variable for the comparison (should you be using k itself). Can you provide valid examples for your starting values?

[Irene Zhou]

I have updated the question. Thank you!

[Andy]

The way the code is written, you set dt =0.08 at the start and from then it only gets smaller. When you add do the sum t(K+1) = dt + t(k) this is tending to value that is less than 11 so the while loop condition will never be met.

Should you be resetting dt=0.08 at the same time as incrementing i ?

[Andy]

It starts going "wrong" on the 33 run through the loop when t=2.3414, I haven't plotted it but what is the function doing at this point?

[Irene Zhou]

At the end of my Else statement, I update k = k+1 so t(k) is able to reach 11 in theory.

The purpose of comparing the Euler's method to the midpoint method is to take dynamic step sizes. If the difference is too big, take a smaller step. If the difference is small enough, try a bigger step, which I haven't got to.

It starts to go wrong when dif hits the threshold I set, 0.1 in this case. dt keeps getting reduced by half while nothing else gets updated, which is what I am trying to figure out.