I get different solutions from ODE45

I tried to solve MIT's open courseware MatLab assigments. On the assigment they wanted plot different colors of lines depending on their value if it bigger than zero it would be red, it isn't it would black.

Here is the differential system and constans values.

My ODE45 function's code.

function dif_system = diferantial(t, x);
% Solver for the differantial system
%   alphan(V), alphah(V), alpham(V), betan(V), betam(V), betah(V) 
%   these are the funciton scripts MIT gave us to use.

n = x(1);
m = x(2);
h = x(3);
V = x(4);

An = alphan(V); Bn = betan(V);
Am = alpham(V); Bm = betam(V);
Ah = alphah(V); Bh = betah(V);

% Membrane Capaticane
C = 1; 

% Conductances
Gk = 36; % Potassium
Gna = 120; % Sodium
Gl = 0.3; % Leak channels

% Reversal potetias 
Ek = -72; % Potassium
Ena = 55; % Sodium
El = -49.4; % Leak channels

dif_system = [(1-n)*An-(n*Bn);
              (1-m)*Am-(m*Bm); 
              (1-h)*Ah-(h*Bh); 
              -C^-1*(Gk*n^4*(V-Ek)+Gna*m^3*h*(V-Ena)+Gl*(V-El))];
   

end

When i run the code my whole values are bigger than the zero(I shoved it down below).

Where i'm doing wong?

[Torsten]

Why are all your values bigger than zero ?

And you did not include the complete code.

[Emir Dönmez]

Firstly, yes i just upload the funciton scripti will give you the whole code down there. Secondly i don't understand the reason why my whole values are bigger than zero.

Qestion

Whole code

%% QUE_5

t_span = [0 20]; % 20ms
n = 0.5; m = 0.5; h = 0.5; V = -60;
conditions = [n m h V];
ySS = [n m h V];

hold on;
title("Approaching Steady State");
xlabel("Time (ms)"), ylabel("Transmembrane Voltage (V)");
for i=1:10
    ySS(4) = ySS(4)+1;
    [t, x] = ode45(@(t, x) diferantial(t, x), t_span, ySS);
    if max(x(:,4)) > 0
        plot(t, x(:,4), "Color", "red", "LineWidth", 1)
    else 
        plot(t, x(:,4), "Color", "black", "LineWidth", 1)
    end
end
hold off

function dif_system = diferantial(t, x);

n = x(1);
m = x(2);
h = x(3);
V = x(4);

An = alphan(V); Bn = betan(V);
Am = alpham(V); Bm = betam(V);
Ah = alphah(V); Bh = betah(V);

C = 1; Gk = 36; Gna = 120; Gl = 0.3;
Ek = -72; Ena = 55; El = -49.4;

dif_system = [(1-n)*An-(n*Bn);
              (1-m)*Am-(m*Bm); 
              (1-h)*Ah-(h*Bh); 
              -C^-1*(Gk*n^4*(V-Ek)+Gna*m^3*h*(V-Ena)+Gl*(V-El))];
   

end
% returns the beta gating parameter for the m gate at voltage V
function b = betam(V)

b = 4*exp(-(V+60)/18);
end
% returns the beta gating parameter for the n gate at voltage V
function b = betan(V)

b = 0.125*exp(-(V+60)/80);
end
% returns the beta gating parameter for the h gate at voltage V
function b = betah(V)

b = 1./(exp(-(V+30)/10)+1);
end
% returns the alpha gating parameter for the h gate at voltage V
function a = alphah(V)

a = 0.07*exp(-0.05*(V+60));
end
% returns the alpha gating parameter for the m gate at voltage V
function a = alpham(V)

badIndices = find(abs(V+35)<1e-4);
goodIndices = find(abs(V+35)>=1e-4);

a(goodIndices) = -0.1*(V(goodIndices)+35)./(exp(-(V(goodIndices)+35)/10)-1);
if(max(size(badIndices))>0)
    a(badIndices) = 1./(1-(V(badIndices)+35)/20);
end
end
% returns the alpha gating parameter for the n gate at voltage V
function a = alphan(V)

badIndices = find(abs(V+50)<1e-4);
goodIndices = find(abs(V+50)>=1e-4);

a(goodIndices) = -0.01*(V(goodIndices)+50)./(exp(-(V(goodIndices)+50)/10)-1);
if(max(size(badIndices))>0)
    a(badIndices) = .1./(1-(V(badIndices)+50)/20);
end
end

I have attached the files that needed with function scritp file.

[Torsten]

With the code you show above, the complete curve x(:,4) will be plotted in red if only one of its values is > 0 (because in this case, the condition max(x(:,4)) > 0 is true).

So reconsider the if-statement.

[Emir Dönmez]

But when i run the code i get 10 lines in graph and all of are bigger than zero? I will get the same plot even if i change the if-statement.

[Torsten]

Yes, the peak value for V is always bigger than 0. Thus the curves are to be plotted in red - as you do in your code.

But according to the assignment, I'd expect that some curves should appear in black. So reconsider your initial conditions and equations.

[Emir Dönmez]

I know some problem in the conditions or equetions that's way i asked.

[Torsten]

Usually, we cannot find modelling errors. We can only help with MATLAB syntax errors and sometimes with solver problems.

Answer 1

[Walter Roberson]

Plotting the whole line in red if the voltage peak for the line is above 0 is what the question asks. It does not ask for the part of the line above 0 to be plotted in red.

Look for example at the wording about zooming in to see the threshold that separates the red lines from the black lines: that is definitely talking about red for the entire lines, not for line segments.

[Emir Dönmez]

I understood on the what they want me on the question, they wanted me to paint depending on peak values but the problem is whole lines peak values are way bigger than the zero.