Need Help Converting Berkeley Madonna code to Matlab code

Here is the Berkeley Madonna code:

{Top model}

{Reservoirs}

d/dt (Q) = + Stimulus - Imemb

INIT Q = -65/cap

{Flows}

Stimulus = Intensity*SquarePulse(3,.5) {at t=3 of 0.5 duration}

Imemb = IL+IK+INa

{Functions}

Intensity = 100 {microamps}

cap = 1

E = Q/cap

{Submodel "INa_"}

{Functions}

ENa = 50

INa = gNa*(E-ENa)

GNaMax = 120

gNa = GNaMax*m*m*m*h

{Submodel "m_gates"}

{Reservoirs}

d/dt (m) = + m_prod - m_decay

INIT m = am/(am+bm)

{Flows}

m_prod = am*(1-m)

m_decay = bm*m

{Functions}

am = 0.1*(E+40)/(1-exp(-(E+40)/10))

bm = 4*exp(-(E+65)/18)

{Submodel "h_gates"}

{Reservoirs}

d/dt (h) = + h_prod - h_decay

INIT h = ah/(ah+bh)

{Flows}

h_prod = ah*(1-h)

h_decay = bh*h

{Functions}

ah = 0.07*exp(-(E+65)/20)

bh = 1/(exp(-(E+35)/10)+1)

{Submodel "IK_"}

{Functions}

EK = -77

IK = gK*(E-EK)

gK_max = 36

gK = gK_max*n*n*n*n

{Submodel "n_gates"}

{Reservoirs}

d/dt (n) = + n_prod - n_decay

INIT n = an/(an+bn)

{Flows}

n_prod = an*(1-n)

n_decay = bn*n

{Functions}

an = 0.01*(E+55)/(1-exp(-(E+55)/10))

bn = 0.125*exp(-(E+65)/80)

{Submodel "IL_"}

{Functions}

IL = gL*(E-EL)

EL = -54.4

gL = .3

{Globals}

{End Globals}