What am I doing wrong while defining this ODE function?

What I'm trying to do:

Write a script called HH.m which will solve this system of ODEs and plot the transmembrane voltage. First, we’ll run the system to steady state. Run the simulation for 20ms (the timescale of the equations is ms) with initial values: n = 0.5; m = 0.5;h = 0.5;V = −60 (ode45). Store the steady-state value of all 4 parameters in a vector called ySS. Make a new figure and plot the timecourse of V (t ) to verify that it reaches steady state by the end of the simulation.

What I did:

Function File:

function [dndt dmdt dhdt dVdt]= HHODEfun(t,V)

C=1; GK=36; GNa=120; GL=0.3; EK=-72; ENa=55; EL=-49.4; %Constants

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

dndt=(1-n).*alphan(Vol)-n.*betan(Vol);

dmdt=(1-m).*alpham(Vol)-m.*betam(Vol);

dhdt=(1-h).*alphah(Vol)-h.*betah(Vol);

dVdt=(-1./C).*((GK.*(n.^4).*(Vol-EK))+GNa.*(m.^3).*h.*(Vol-ENa)+GL.*(Vol-EL));

Script:

clc

[t n m h V]=ode45(@HHODEfun,[0 0.02],[0.5,0.5,0.5,-60]);

ySS=[n m h V];

figure

plot(t,V)

which is obviously not correct. what am I doing wrong?

Answer 1

[Star Strider]

Assuming the rest of your code is correct (I didn’t check), you need to change it to return the vector of derivatives as a single column vector, not individual derivatives:

function derivs = HHODEfun(t,V)

C=1; GK=36; GNa=120; GL=0.3; EK=-72; ENa=55; EL=-49.4; %Constants

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

dndt=(1-n).*alphan(Vol)-n.*betan(Vol);
dmdt=(1-m).*alpham(Vol)-m.*betam(Vol);
dhdt=(1-h).*alphah(Vol)-h.*betah(Vol);
dVdt=(-1./C).*((GK.*(n.^4).*(Vol-EK))+GNa.*(m.^3).*h.*(Vol-ENa)+GL.*(Vol-EL));
derivs = [dndt dmdt dhdt dVdt]';
end

and your ode45 call becomes:

[t, V]=ode45(@HHODEfun,[0 0.02],[0.5,0.5,0.5,-60]);

NOTE: I could not run this to test it because alphan and the others are not defined. Define them and it should work.