How to resolve this problem (differential equations) ?

edit

Hi ! I wrote a code for a project which concerns differential equations but the problem is I don't have the good curves, even if the code is right (I think). Or it's a problem of syntax but I spend a lot of time searching this, without success. I give you the code accompanied with the formulas and the exit of the code. I can't add hyperlink so I write like this or you can't understand my problem. Sorry for inconvenience !

http://image.noelshack.com/fichiers/2...

http://image.noelshack.com/fichiers/2...

http://image.noelshack.com/fichiers/2...

I'm supposed to get them curves :/ The only problem seems to be black curve (on SageMath) which is the green curve on Excel (third picture) and I don't wanted to put the black curve of Excel (so it's normal if it doesn't appear).

http://image.noelshack.com/fichiers/2...

# y[0] = X (Number of likely ie those not infected) = B - nu*X - c*lmd*X
# y[1] = Y (Number of contagious) = c*X*lmd - Y*(v + nu)
# y[2] = Z (Number of non-infectious HIV) = v*Y*(1 - p) - nu*Z
# y[3] = A (Number of AIDS patients) = p*v*Y - A*(d + nu)
# y[4] = lmd (Proportion of contagious, which made themselves pass HIV on the total population) = (Y*beta)/(X + Y + Z)

# params[0] = nu (Natural death rate) == 0.03125 car nu + d ~= d (1 à 1.33/an)
# params[1] = d (Death rates by disease) == 1 (1/an)
# params[2] = B (Nombre de susceptibles immigrants) == 13333.3/an
# params[3] = v (Number of immigrants may) == 0.2/an
# params[4] = c (Number of sexual partners) == 48 (2 à 6/mois soit 48/an)
# params[5] = p (Proportion of catching HIV) == 0.3 (10 à 30%)
# params[6] = beta (Transmission probability) == 0.0014
# t (Time)

@interact
def interactive_function ( nu = slider (0.001, 0.1, 0.005, default = 0.03125),
                           d = slider (0, 1.5, 0.05, default = 1),
                           B = slider (0, 30000, 500, default = 13333.3),
                           v = slider (0, 1, 0.05, default = 0.2),
                           c = slider (0, 365, 1, default = 48),
                           p = slider (0, 1, 0.05, default = 0.3),
                           beta = slider (0, 1, 0.02, default = 0.014) ) :

  def f_1 (t,y,params) :

    return [ params[2] - params[0]*y[0] - params[4]*y[4]*y[0],
             params[4]*y[0]*y[4] - y[1]*(params[3] + params[0]),
             params[3]*y[1]*(1-params[5]) - params[0]*y[2],
             params[5]*params[3]*y[1] - y[3]*(params[1] + params[0]),
             (y[1]*params[6])/(y[0] + y[1] + y[2]) ]


  T = ode_solver()
  T.function = f_1
  T.algorithm = "rk8pd"
  T.ode_solve (y_0 = [60000,40000,0,0,0.0056],
               t_span = [0,35],
               params = [nu,d,B,v,c,p,beta],
               num_points = 1000)

  f = T.solution
  X = [(x[0],x[1][0]) for x in f]
  Y = [(x[0],x[1][1]) for x in f]
  Z = [(x[0],x[1][2]) for x ...