Having problems reproducing Matlab results in Python for a coupled system of ODE's
Question: I just started to use Python and have a bunch of code I want to transfer over
from Matlab. I started with a simple coupled diff eq and can't seem to figure
out what I'm doing wrong. It appears that the second diff eq is almost giving
me the appropriate results but the first equation is just remaining at its
initial condition for the entire time span. The first image is the Python
result and the second is the correct result from Matlab.

1) ![](http://i.imgur.com/VSFpCT3.png) 2) ![](http://i.imgur.com/w7sfWko.png)

    
    
    from scipy.integrate import odeint
    from pylab import *
    
    minfH_ab = lambda v: 1/(1 + exp((v+70)/6))
    taumH_ab = lambda v: 272 + 1499/(1 + exp(-(v+42.2)/8.73))
    
    
    Csn = 9    
    I_ab_sn  = 0
    gL_ab_sn = 0.045 
    El_ab_sn = -50
    gH_ab   = 0.054
    Eh_ab   = -20
    
    def dy_dt(y, t):
        dy1 = (1/Csn)*(I_ab_sn -((gH_ab*y[1]*(y[0]-Eh_ab))+(gL_ab_sn*(y[0]-El_ab_sn))))
        dy2 = (minfH_ab(y[0])-y[1])/taumH_ab(y[0])
        return [dy1, dy2]
    
    
    t = linspace(0,1000,10000)
    y_init = [-50, .0004]
    
    sol = odeint(dy_dt, y_init, t)
    S0 = sol[:, 0]
    S1 = sol[:, 1]
    
    figure()
    plot(t, S0)
    xlabel('time')
    ylabel('voltage')
    title('H & L Current')
    


Answer: Since this turned out to be the answer:

# Beware integer arithmetic.

If you're porting from a language that does not do this by default, you may be
surprised when

    
    
    x=50
    y=6
    z=5
    
    x*y/z == 60
    x*(y/z) == 50
    x/y*z == 40
    x*z/y == 41
    

And so on. Try changing all constants to doubles (`5` -> `5.0`, etc.), and see
if that helps.

EDIT: @BenDundee pointed out that python will be switching syntax to `/`
meaning "real division", and integer division being specified with the `//`
operator. You can switch to that behaviour now, with the line `from __future__
import division`, which should also resolve the problem.

