# Size of variable arrays:
sizeAlgebraic = 82
sizeStates = 29
sizeConstants = 63
from math import *
from numpy import *

def createLegends():
    legend_states = [""] * sizeStates
    legend_rates = [""] * sizeStates
    legend_algebraic = [""] * sizeAlgebraic
    legend_voi = ""
    legend_constants = [""] * sizeConstants
    legend_voi = "time in component environment (second)"
    legend_states[0] = "V in component membrane (millivolt)"
    legend_constants[0] = "R in component membrane (joule_per_kilomole_kelvin)"
    legend_constants[1] = "T in component membrane (kelvin)"
    legend_constants[2] = "F in component membrane (coulomb_per_mole)"
    legend_constants[3] = "C in component membrane (microF)"
    legend_constants[53] = "RTONF in component membrane (millivolt)"
    legend_algebraic[66] = "i_Na in component fast_sodium_current (nanoA)"
    legend_algebraic[71] = "i_CaL in component L_type_calcium_current (nanoA)"
    legend_algebraic[67] = "i_to in component transient_outward_potassium_current (nanoA)"
    legend_algebraic[42] = "i_Kr in component rapid_delayed_rectifier_potassium_current (nanoA)"
    legend_algebraic[34] = "i_f in component hyperpolarising_activated_current (nanoA)"
    legend_algebraic[68] = "i_st in component sustained_outward_potassium_current (nanoA)"
    legend_algebraic[50] = "i_K1 in component time_independent_potassium_current (nanoA)"
    legend_algebraic[65] = "i_NaCa in component sodium_calcium_exchange_current (nanoA)"
    legend_algebraic[52] = "i_p in component sodium_potassium_pump (nanoA)"
    legend_algebraic[51] = "i_b in component background_current (nanoA)"
    legend_algebraic[70] = "i_ACh in component acetylcholine_sensitive_current (nanoA)"
    legend_algebraic[17] = "i_Stim in component membrane (nanoA)"
    legend_constants[4] = "stim_start in component membrane (second)"
    legend_constants[5] = "stim_end in component membrane (second)"
    legend_constants[6] = "stim_period in component membrane (second)"
    legend_constants[7] = "stim_duration in component membrane (second)"
    legend_constants[8] = "stim_amplitude in component membrane (nanoA)"
    legend_constants[9] = "g_f in component hyperpolarising_activated_current (microS)"
    legend_constants[10] = "ACh in component acetylcholine_sensitive_current (millimolar)"
    legend_states[1] = "y in component hyperpolarising_activated_current_y_gate (dimensionless)"
    legend_algebraic[0] = "y_inf in component hyperpolarising_activated_current_y_gate (dimensionless)"
    legend_algebraic[19] = "tau_y in component hyperpolarising_activated_current_y_gate (second)"
    legend_constants[11] = "g_Kr in component rapid_delayed_rectifier_potassium_current (microS)"
    legend_constants[55] = "E_K in component rapid_delayed_rectifier_potassium_current (millivolt)"
    legend_constants[12] = "Ki in component intracellular_potassium_concentration (millimolar)"
    legend_constants[13] = "Kc in component extracellular_potassium_concentration (millimolar)"
    legend_states[2] = "paf in component rapid_delayed_rectifier_potassium_current_paf_gate (dimensionless)"
    legend_states[3] = "pas in component rapid_delayed_rectifier_potassium_current_pas_gate (dimensionless)"
    legend_states[4] = "pik in component rapid_delayed_rectifier_potassium_current_pik_gate (dimensionless)"
    legend_algebraic[1] = "paf_infinity in component rapid_delayed_rectifier_potassium_current_paf_gate (dimensionless)"
    legend_algebraic[20] = "tau_paf in component rapid_delayed_rectifier_potassium_current_paf_gate (second)"
    legend_algebraic[2] = "pas_infinity in component rapid_delayed_rectifier_potassium_current_pas_gate (dimensionless)"
    legend_algebraic[21] = "tau_pas in component rapid_delayed_rectifier_potassium_current_pas_gate (second)"
    legend_algebraic[3] = "pik_infinity in component rapid_delayed_rectifier_potassium_current_pik_gate (dimensionless)"
    legend_algebraic[22] = "alpha_pik in component rapid_delayed_rectifier_potassium_current_pik_gate (per_second)"
    legend_algebraic[35] = "beta_pik in component rapid_delayed_rectifier_potassium_current_pik_gate (per_second)"
    legend_algebraic[43] = "tau_pik in component rapid_delayed_rectifier_potassium_current_pik_gate (second)"
    legend_constants[14] = "g_K1 in component time_independent_potassium_current (microS)"
    legend_algebraic[49] = "g_K1_prime in component time_independent_potassium_current (microS)"
    legend_constants[15] = "g_b in component background_current (microS)"
    legend_constants[16] = "E_b in component background_current (millivolt)"
    legend_constants[17] = "I_p in component sodium_potassium_pump (nanoA)"
    legend_constants[18] = "Nai in component intracellular_sodium_concentration (millimolar)"
    legend_constants[19] = "kNaCa in component sodium_calcium_exchange_current (nanoA)"
    legend_algebraic[61] = "x1 in component sodium_calcium_exchange_current (dimensionless)"
    legend_algebraic[62] = "x2 in component sodium_calcium_exchange_current (dimensionless)"
    legend_algebraic[63] = "x3 in component sodium_calcium_exchange_current (dimensionless)"
    legend_algebraic[64] = "x4 in component sodium_calcium_exchange_current (dimensionless)"
    legend_algebraic[58] = "k41 in component sodium_calcium_exchange_current (dimensionless)"
    legend_constants[54] = "k34 in component sodium_calcium_exchange_current (dimensionless)"
    legend_algebraic[55] = "k23 in component sodium_calcium_exchange_current (dimensionless)"
    legend_algebraic[56] = "k21 in component sodium_calcium_exchange_current (dimensionless)"
    legend_algebraic[54] = "k32 in component sodium_calcium_exchange_current (dimensionless)"
    legend_constants[58] = "k43 in component sodium_calcium_exchange_current (dimensionless)"
    legend_algebraic[60] = "k12 in component sodium_calcium_exchange_current (dimensionless)"
    legend_algebraic[59] = "k14 in component sodium_calcium_exchange_current (dimensionless)"
    legend_constants[20] = "Qci in component sodium_calcium_exchange_current (dimensionless)"
    legend_constants[21] = "Qn in component sodium_calcium_exchange_current (dimensionless)"
    legend_constants[22] = "Qco in component sodium_calcium_exchange_current (dimensionless)"
    legend_constants[23] = "Kci in component sodium_calcium_exchange_current (millimolar)"
    legend_constants[24] = "K1ni in component sodium_calcium_exchange_current (millimolar)"
    legend_constants[25] = "K2ni in component sodium_calcium_exchange_current (millimolar)"
    legend_constants[26] = "K3ni in component sodium_calcium_exchange_current (millimolar)"
    legend_constants[27] = "Kcni in component sodium_calcium_exchange_current (millimolar)"
    legend_constants[28] = "K3no in component sodium_calcium_exchange_current (millimolar)"
    legend_constants[29] = "K1no in component sodium_calcium_exchange_current (millimolar)"
    legend_constants[30] = "K2no in component sodium_calcium_exchange_current (millimolar)"
    legend_constants[31] = "Kco in component sodium_calcium_exchange_current (millimolar)"
    legend_algebraic[53] = "do in component sodium_calcium_exchange_current (dimensionless)"
    legend_algebraic[57] = "di in component sodium_calcium_exchange_current (dimensionless)"
    legend_constants[32] = "Cao in component extracellular_calcium_concentration (millimolar)"
    legend_constants[33] = "Nao in component extracellular_sodium_concentration (millimolar)"
    legend_states[5] = "Casub in component intracellular_calcium_concentration (millimolar)"
    legend_constants[34] = "g_Na in component fast_sodium_current (microlitre_per_second)"
    legend_constants[56] = "E_Na in component fast_sodium_current (millivolt)"
    legend_states[6] = "m in component fast_sodium_current_m_gate (dimensionless)"
    legend_states[7] = "h1 in component fast_sodium_current_h1_gate (dimensionless)"
    legend_states[8] = "h2 in component fast_sodium_current_h2_gate (dimensionless)"
    legend_algebraic[23] = "alpha_m in component fast_sodium_current_m_gate (per_second)"
    legend_algebraic[36] = "beta_m in component fast_sodium_current_m_gate (per_second)"
    legend_constants[35] = "delta_m in component fast_sodium_current_m_gate (millivolt)"
    legend_algebraic[4] = "E0_m in component fast_sodium_current_m_gate (millivolt)"
    legend_algebraic[5] = "alpha_h1 in component fast_sodium_current_h1_gate (per_second)"
    legend_algebraic[24] = "beta_h1 in component fast_sodium_current_h1_gate (per_second)"
    legend_algebraic[37] = "h1_inf in component fast_sodium_current_h1_gate (dimensionless)"
    legend_algebraic[44] = "tau_h1 in component fast_sodium_current_h1_gate (second)"
    legend_algebraic[6] = "alpha_h2 in component fast_sodium_current_h2_gate (per_second)"
    legend_algebraic[25] = "beta_h2 in component fast_sodium_current_h2_gate (per_second)"
    legend_algebraic[38] = "h2_inf in component fast_sodium_current_h2_gate (dimensionless)"
    legend_algebraic[45] = "tau_h2 in component fast_sodium_current_h2_gate (second)"
    legend_constants[36] = "g_CaL in component L_type_calcium_current (microS)"
    legend_constants[37] = "E_CaL in component L_type_calcium_current (millivolt)"
    legend_states[9] = "d in component L_type_calcium_current_d_gate (dimensionless)"
    legend_states[10] = "f in component L_type_calcium_current_f_gate (dimensionless)"
    legend_states[11] = "f2 in component L_type_calcium_current_f2_gate (dimensionless)"
    legend_algebraic[7] = "alpha_d in component L_type_calcium_current_d_gate (per_second)"
    legend_algebraic[26] = "beta_d in component L_type_calcium_current_d_gate (per_second)"
    legend_algebraic[39] = "d_inf in component L_type_calcium_current_d_gate (dimensionless)"
    legend_algebraic[46] = "tau_d in component L_type_calcium_current_d_gate (second)"
    legend_constants[38] = "act_shift in component L_type_calcium_current_d_gate (millivolt)"
    legend_constants[39] = "slope_factor_act in component L_type_calcium_current_d_gate (millivolt)"
    legend_algebraic[8] = "f_inf in component L_type_calcium_current_f_gate (dimensionless)"
    legend_algebraic[27] = "tau_f in component L_type_calcium_current_f_gate (second)"
    legend_constants[40] = "inact_shift in component L_type_calcium_current_f_gate (millivolt)"
    legend_algebraic[9] = "f2_inf in component L_type_calcium_current_f2_gate (dimensionless)"
    legend_algebraic[28] = "tau_f2 in component L_type_calcium_current_f2_gate (second)"
    legend_constants[41] = "inact_shift in component L_type_calcium_current_f2_gate (millivolt)"
    legend_constants[57] = "E_K in component transient_outward_potassium_current (millivolt)"
    legend_constants[42] = "g_to in component transient_outward_potassium_current (microS)"
    legend_states[12] = "r in component transient_outward_potassium_current_r_gate (dimensionless)"
    legend_states[13] = "q_fast in component transient_outward_potassium_current_qfast_gate (dimensionless)"
    legend_states[14] = "q_slow in component transient_outward_potassium_current_qslow_gate (dimensionless)"
    legend_algebraic[29] = "tau_r in component transient_outward_potassium_current_r_gate (second)"
    legend_algebraic[10] = "r_infinity in component transient_outward_potassium_current_r_gate (dimensionless)"
    legend_algebraic[30] = "tau_qfast in component transient_outward_potassium_current_qfast_gate (second)"
    legend_algebraic[11] = "qfast_infinity in component transient_outward_potassium_current_qfast_gate (dimensionless)"
    legend_algebraic[31] = "tau_qslow in component transient_outward_potassium_current_qslow_gate (second)"
    legend_algebraic[12] = "qslow_infinity in component transient_outward_potassium_current_qslow_gate (dimensionless)"
    legend_constants[43] = "E_st in component sustained_outward_potassium_current (millivolt)"
    legend_constants[44] = "g_st in component sustained_outward_potassium_current (microS)"
    legend_states[15] = "qa in component sustained_outward_potassium_current_qa_gate (dimensionless)"
    legend_states[16] = "qi in component sustained_outward_potassium_current_qi_gate (dimensionless)"
    legend_algebraic[47] = "tau_qa in component sustained_outward_potassium_current_qa_gate (second)"
    legend_algebraic[13] = "qa_infinity in component sustained_outward_potassium_current_qa_gate (dimensionless)"
    legend_algebraic[32] = "alpha_qa in component sustained_outward_potassium_current_qa_gate (per_second)"
    legend_algebraic[40] = "beta_qa in component sustained_outward_potassium_current_qa_gate (per_second)"
    legend_algebraic[48] = "tau_qi in component sustained_outward_potassium_current_qi_gate (second)"
    legend_algebraic[14] = "alpha_qi in component sustained_outward_potassium_current_qi_gate (per_second)"
    legend_algebraic[33] = "beta_qi in component sustained_outward_potassium_current_qi_gate (per_second)"
    legend_algebraic[41] = "qi_infinity in component sustained_outward_potassium_current_qi_gate (dimensionless)"
    legend_algebraic[69] = "g_ACh in component acetylcholine_sensitive_current (microS)"
    legend_constants[45] = "g_ACh_max in component acetylcholine_sensitive_current (microS)"
    legend_constants[46] = "K_ACh in component acetylcholine_sensitive_current (millimolar)"
    legend_states[17] = "achf in component acetylcholine_sensitive_current_achf_gate (dimensionless)"
    legend_states[18] = "achs in component acetylcholine_sensitive_current_achs_gate (dimensionless)"
    legend_constants[47] = "alpha_achf in component acetylcholine_sensitive_current_achf_gate (per_second)"
    legend_algebraic[15] = "beta_achf in component acetylcholine_sensitive_current_achf_gate (per_second)"
    legend_constants[48] = "alpha_achs in component acetylcholine_sensitive_current_achs_gate (per_second)"
    legend_algebraic[16] = "beta_achs in component acetylcholine_sensitive_current_achs_gate (per_second)"
    legend_states[19] = "Cai in component intracellular_calcium_concentration (millimolar)"
    legend_constants[59] = "V_up in component intracellular_calcium_concentration (micrometre3)"
    legend_constants[60] = "V_rel in component intracellular_calcium_concentration (micrometre3)"
    legend_constants[61] = "V_sub in component intracellular_calcium_concentration (micrometre3)"
    legend_constants[62] = "Vi in component intracellular_calcium_concentration (micrometre3)"
    legend_constants[49] = "V_cell in component intracellular_calcium_concentration (micrometre3)"
    legend_algebraic[73] = "i_up in component intracellular_calcium_concentration (millimolar_per_second)"
    legend_algebraic[74] = "i_tr in component intracellular_calcium_concentration (millimolar_per_second)"
    legend_algebraic[76] = "i_rel in component intracellular_calcium_concentration (millimolar_per_second)"
    legend_algebraic[72] = "i_diff in component intracellular_calcium_concentration (millimolar_per_second)"
    legend_states[20] = "Ca_up in component intracellular_calcium_concentration (millimolar)"
    legend_states[21] = "Ca_rel in component intracellular_calcium_concentration (millimolar)"
    legend_constants[50] = "P_rel in component intracellular_calcium_concentration (per_second)"
    legend_constants[51] = "K_up in component intracellular_calcium_concentration (millimolar)"
    legend_constants[52] = "tau_tr in component intracellular_calcium_concentration (second)"
    legend_states[22] = "f_TC in component intracellular_calcium_concentration (dimensionless)"
    legend_states[23] = "f_TMC in component intracellular_calcium_concentration (dimensionless)"
    legend_states[24] = "f_TMM in component intracellular_calcium_concentration (dimensionless)"
    legend_states[25] = "f_CMi in component intracellular_calcium_concentration (dimensionless)"
    legend_states[26] = "f_CMs in component intracellular_calcium_concentration (dimensionless)"
    legend_states[27] = "f_CQ in component intracellular_calcium_concentration (dimensionless)"
    legend_states[28] = "f_CSL in component intracellular_calcium_concentration (dimensionless)"
    legend_algebraic[75] = "diff_f_TC in component intracellular_calcium_concentration (per_second)"
    legend_algebraic[77] = "diff_f_TMC in component intracellular_calcium_concentration (per_second)"
    legend_algebraic[18] = "diff_f_TMM in component intracellular_calcium_concentration (per_second)"
    legend_algebraic[78] = "diff_f_CMi in component intracellular_calcium_concentration (per_second)"
    legend_algebraic[79] = "diff_f_CMs in component intracellular_calcium_concentration (per_second)"
    legend_algebraic[80] = "diff_f_CQ in component intracellular_calcium_concentration (per_second)"
    legend_algebraic[81] = "diff_f_CSL in component intracellular_calcium_concentration (per_second)"
    legend_rates[0] = "d/dt V in component membrane (millivolt)"
    legend_rates[1] = "d/dt y in component hyperpolarising_activated_current_y_gate (dimensionless)"
    legend_rates[2] = "d/dt paf in component rapid_delayed_rectifier_potassium_current_paf_gate (dimensionless)"
    legend_rates[3] = "d/dt pas in component rapid_delayed_rectifier_potassium_current_pas_gate (dimensionless)"
    legend_rates[4] = "d/dt pik in component rapid_delayed_rectifier_potassium_current_pik_gate (dimensionless)"
    legend_rates[6] = "d/dt m in component fast_sodium_current_m_gate (dimensionless)"
    legend_rates[7] = "d/dt h1 in component fast_sodium_current_h1_gate (dimensionless)"
    legend_rates[8] = "d/dt h2 in component fast_sodium_current_h2_gate (dimensionless)"
    legend_rates[9] = "d/dt d in component L_type_calcium_current_d_gate (dimensionless)"
    legend_rates[10] = "d/dt f in component L_type_calcium_current_f_gate (dimensionless)"
    legend_rates[11] = "d/dt f2 in component L_type_calcium_current_f2_gate (dimensionless)"
    legend_rates[12] = "d/dt r in component transient_outward_potassium_current_r_gate (dimensionless)"
    legend_rates[13] = "d/dt q_fast in component transient_outward_potassium_current_qfast_gate (dimensionless)"
    legend_rates[14] = "d/dt q_slow in component transient_outward_potassium_current_qslow_gate (dimensionless)"
    legend_rates[15] = "d/dt qa in component sustained_outward_potassium_current_qa_gate (dimensionless)"
    legend_rates[16] = "d/dt qi in component sustained_outward_potassium_current_qi_gate (dimensionless)"
    legend_rates[17] = "d/dt achf in component acetylcholine_sensitive_current_achf_gate (dimensionless)"
    legend_rates[18] = "d/dt achs in component acetylcholine_sensitive_current_achs_gate (dimensionless)"
    legend_rates[20] = "d/dt Ca_up in component intracellular_calcium_concentration (millimolar)"
    legend_rates[21] = "d/dt Ca_rel in component intracellular_calcium_concentration (millimolar)"
    legend_rates[19] = "d/dt Cai in component intracellular_calcium_concentration (millimolar)"
    legend_rates[5] = "d/dt Casub in component intracellular_calcium_concentration (millimolar)"
    legend_rates[22] = "d/dt f_TC in component intracellular_calcium_concentration (dimensionless)"
    legend_rates[23] = "d/dt f_TMC in component intracellular_calcium_concentration (dimensionless)"
    legend_rates[24] = "d/dt f_TMM in component intracellular_calcium_concentration (dimensionless)"
    legend_rates[25] = "d/dt f_CMi in component intracellular_calcium_concentration (dimensionless)"
    legend_rates[26] = "d/dt f_CMs in component intracellular_calcium_concentration (dimensionless)"
    legend_rates[27] = "d/dt f_CQ in component intracellular_calcium_concentration (dimensionless)"
    legend_rates[28] = "d/dt f_CSL in component intracellular_calcium_concentration (dimensionless)"
    return (legend_states, legend_algebraic, legend_voi, legend_constants)

def initConsts():
    constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
    states[0] = -69.760276376489
    constants[0] = 8314.472
    constants[1] = 310
    constants[2] = 96485.3415
    constants[3] = 2.9e-5
    constants[4] = 0.1
    constants[5] = 99999
    constants[6] = 1
    constants[7] = 0.001
    constants[8] = -2
    constants[9] = 0
    constants[10] = 0
    states[1] = 0.19584111039096
    constants[11] = 0.002
    constants[12] = 140
    constants[13] = 5.4
    states[2] = 0.00141762995766447
    states[3] = 0.00539771950846456
    states[4] = 0.98638204514681
    constants[14] = 0.015
    constants[15] = 0.002
    constants[16] = -40
    constants[17] = 0.1968
    constants[18] = 8
    constants[19] = 5.916
    constants[20] = 0.1369
    constants[21] = 0.4315
    constants[22] = 0
    constants[23] = 0.0207
    constants[24] = 395.3
    constants[25] = 2.289
    constants[26] = 26.44
    constants[27] = 26.44
    constants[28] = 4.663
    constants[29] = 1628
    constants[30] = 561.4
    constants[31] = 3.663
    constants[32] = 2
    constants[33] = 140
    states[5] = 3.27335718697622e-5
    constants[34] = 5e-7
    states[6] = 0.0132200747771872
    states[7] = 0.706622937059237
    states[8] = 0.701626826712569
    constants[35] = 1e-5
    constants[36] = 0.021
    constants[37] = 62.1
    states[9] = 4.23500474189711e-5
    states[10] = 0.998435073735753
    states[11] = 0.998424216938754
    constants[38] = 0
    constants[39] = -6.61
    constants[40] = -5
    constants[41] = -5
    constants[42] = 0.014
    states[12] = 0.00894826428663828
    states[13] = 0.994837524153424
    states[14] = 0.427382372349565
    constants[43] = -37.4
    constants[44] = 0
    states[15] = 0.0910882041816457
    states[16] = 0.890389014175329
    constants[45] = 0.0198
    constants[46] = 0.00035
    states[17] = 0.74242774587522
    states[18] = 0.746918323597392
    constants[47] = 73.1
    constants[48] = 3.7
    states[19] = 3.73317163732586e-5
    constants[49] = 4.39823e-6
    states[20] = 0.818671555213184
    states[21] = 0.682159360306652
    constants[50] = 1805.6
    constants[51] = 0.0006
    constants[52] = 0.06
    states[22] = 0.00739583869345967
    states[23] = 0.133773505989393
    states[24] = 0.765246381247448
    states[25] = 0.0154714370092264
    states[26] = 0.0135767851016544
    states[27] = 0.449992033196647
    states[28] = 1.21722016147587e-5
    constants[53] = (constants[0]*constants[1])/constants[2]
    constants[54] = constants[33]/(constants[28]+constants[33])
    constants[55] = constants[53]*log(constants[13]/constants[12])
    constants[56] = constants[53]*log(constants[33]/constants[18])
    constants[57] = constants[53]*log(constants[13]/constants[12])
    constants[58] = constants[18]/(constants[26]+constants[18])
    constants[59] = 0.0116000*constants[49]
    constants[60] = 0.00120000*constants[49]
    constants[61] = 0.0100000*constants[49]
    constants[62] = 0.460000*constants[49]-constants[61]
    return (states, constants)

def computeRates(voi, states, constants):
    rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic
    algebraic[15] = 120.000/(1.00000+exp(-(states[0]+50.0000)/15.0000))
    rates[17] = constants[47]*(1.00000-states[17])-algebraic[15]*states[17]
    algebraic[16] = 5.82000/(1.00000+exp(-(states[0]+50.0000)/15.0000))
    rates[18] = constants[48]*(1.00000-states[18])-algebraic[16]*states[18]
    algebraic[18] = 2277.00*2.50000*((1.00000-states[23])-states[24])-751.000*states[24]
    rates[24] = algebraic[18]
    algebraic[0] = 1.00000/(1.00000+exp(((states[0]+83.1900)-(-7.20000*(power(constants[10], 0.690000)))/(power(1.26000e-05, 0.690000)+power(constants[10], 0.690000)))/13.5600))
    algebraic[19] = 0.250000+2.00000*exp(-(power(states[0]+70.0000, 2.00000))/500.000)
    rates[1] = (algebraic[0]-states[1])/algebraic[19]
    algebraic[1] = 1.00000/(1.00000+exp((states[0]+10.2200)/-8.50000))
    algebraic[20] = 1.00000/(17.0000*exp(0.0398000*states[0])+0.211000*exp(-0.0510000*states[0]))
    rates[2] = (algebraic[1]-states[2])/algebraic[20]
    algebraic[2] = 1.00000/(1.00000+exp((states[0]+10.2200)/-8.50000))
    algebraic[21] = 0.335810+0.906730*exp(-(power(states[0]+10.0000, 2.00000))/988.050)
    rates[3] = (algebraic[2]-states[3])/algebraic[21]
    algebraic[8] = 1.00000/(1.00000+exp((states[0]-(-24.0000+constants[40]))/6.31000))
    algebraic[27] = 0.0100000+0.153900*exp(-(power(states[0]+40.0000, 2.00000))/185.670)
    rates[10] = (algebraic[8]-states[10])/algebraic[27]
    algebraic[9] = 1.00000/(1.00000+exp((states[0]-(-24.0000+constants[41]))/6.31000))
    algebraic[28] = 0.0600000+0.480760*2.25000*exp(-(power(states[0]--40.0000, 2.00000))/138.040)
    rates[11] = (algebraic[9]-states[11])/algebraic[28]
    algebraic[29] = 0.000596000+0.00311800/(1.03700*exp(0.0900000*(states[0]+30.6100))+0.396000*exp(-0.120000*(states[0]+23.8400)))
    algebraic[10] = 1.00000/(1.00000+exp((states[0]-7.44000)/-16.4000))
    rates[12] = (algebraic[10]-states[12])/algebraic[29]
    algebraic[30] = 0.0126600+4.72716/(1.00000+exp((states[0]+154.500)/23.9600))
    algebraic[11] = 1.00000/(1.00000+exp((states[0]+33.8000)/6.12000))
    rates[13] = (algebraic[11]-states[13])/algebraic[30]
    algebraic[31] = 0.100000+4.00000*exp(-(power(states[0]+65.0000, 2.00000))/500.000)
    algebraic[12] = 1.00000/(1.00000+exp((states[0]+33.8000)/6.12000))
    rates[14] = (algebraic[12]-states[14])/algebraic[31]
    algebraic[4] = states[0]+44.4000
    algebraic[23] = custom_piecewise([less(fabs(algebraic[4]) , constants[35]), (-460.000*-12.6730)/exp(algebraic[4]/-12.6730) , True, (-460.000*algebraic[4])/(exp(algebraic[4]/-12.6730)-1.00000)])
    algebraic[36] = 18400.0*exp(algebraic[4]/-12.6730)
    rates[6] = algebraic[23]*(1.00000-states[6])-algebraic[36]*states[6]
    algebraic[3] = (1.00000/(1.00000+exp((states[0]+4.90000)/15.1400)))*(1.00000-0.300000*exp(-(power(states[0], 2.00000))/500.000))
    algebraic[22] = 92.0100*exp(-0.0183000*states[0])
    algebraic[35] = 603.600*exp(0.00942000*states[0])
    algebraic[43] = 1.00000/(algebraic[22]+algebraic[35])
    rates[4] = (algebraic[3]-states[4])/algebraic[43]
    algebraic[5] = 44.9000*exp((states[0]+66.9000)/-5.57000)
    algebraic[24] = 1491.00/(1.00000+323.300*exp((states[0]+94.6000)/-12.9000))
    algebraic[37] = algebraic[5]/(algebraic[5]+algebraic[24])
    algebraic[44] = 0.0300000/(1.00000+exp((states[0]+40.0000)/6.00000))+0.000350000
    rates[7] = (algebraic[37]-states[7])/algebraic[44]
    algebraic[6] = 44.9000*exp((states[0]+66.9000)/-5.57000)
    algebraic[25] = 1491.00/(1.00000+323.300*exp((states[0]+94.6000)/-12.9000))
    algebraic[38] = algebraic[6]/(algebraic[6]+algebraic[25])
    algebraic[45] = 0.120000/(1.00000+exp((states[0]+60.0000)/2.00000))+0.00295000
    rates[8] = (algebraic[38]-states[8])/algebraic[45]
    algebraic[39] = 1.00000/(1.00000+exp((states[0]-(-3.20000+constants[38]))/constants[39]))
    algebraic[7] = (-26.1200*(states[0]+35.0000))/(exp((states[0]+35.0000)/-2.50000)-1.00000)+(-78.1100*states[0])/(exp(-0.208000*states[0])-1.00000)
    algebraic[26] = (10.5200*(states[0]-5.00000))/(exp(0.400000*(states[0]-5.00000))-1.00000)
    algebraic[46] = 1.00000/(algebraic[7]+algebraic[26])
    rates[9] = (algebraic[39]-states[9])/algebraic[46]
    algebraic[32] = 1.00000/(0.150000*exp(-states[0]/11.0000)+0.200000*exp(-states[0]/700.000))
    algebraic[40] = 1.00000/(16.0000*exp(states[0]/8.00000)+15.0000*exp(states[0]/50.0000))
    algebraic[47] = 0.00100000/(algebraic[32]+algebraic[40])
    algebraic[13] = 1.00000/(1.00000+exp((states[0]--49.1000)/-8.98000))
    rates[15] = (algebraic[13]-states[15])/algebraic[47]
    algebraic[14] = 0.150400/(3100.00*exp(states[0]/13.0000)+700.000*exp(states[0]/70.0000))
    algebraic[33] = 0.150400/(95.0000*exp(-states[0]/10.0000)+50.0000*exp(-states[0]/700.000))+0.000229000/(1.00000+exp(-states[0]/5.00000))
    algebraic[48] = 0.00100000/(algebraic[14]+algebraic[33])
    algebraic[41] = algebraic[14]/(algebraic[14]+algebraic[33])
    rates[16] = (algebraic[41]-states[16])/algebraic[48]
    algebraic[66] = (((constants[34]*(power(states[6], 3.00000))*(0.635000*states[7]+0.365000*states[8])*constants[33]*states[0]*constants[2])/constants[53])*(exp((states[0]-constants[56])/constants[53])-1.00000))/(exp(states[0]/constants[53])-1.00000)
    algebraic[69] = (constants[45]*states[17]*states[18]*(power(constants[10], 1.50000)))/(power(constants[46], 1.50000)+power(constants[10], 1.50000))
    algebraic[70] = (((algebraic[69]*constants[13])/(10.0000+constants[13]))*(states[0]-constants[55]))/(1.00000+exp(((states[0]-constants[55])-140.000)/(2.50000*constants[53])))
    algebraic[71] = constants[36]*states[9]*(0.675000*states[10]+0.325000*states[11])*(states[0]-constants[37])*(1.00000-((algebraic[70]*constants[10])/(9.00000e-05+constants[10]))/1.00000)
    algebraic[67] = constants[42]*states[12]*(0.450000*states[13]+0.550000*states[14])*(states[0]-constants[57])
    algebraic[42] = constants[11]*(0.900000*states[2]+0.100000*states[3])*states[4]*(states[0]-constants[55])
    algebraic[34] = states[1]*constants[9]*(states[0]--30.0000)
    algebraic[68] = constants[44]*states[15]*states[16]*(states[0]-constants[43])
    algebraic[49] = constants[14]*(0.500000+0.500000/(1.00000+exp((states[0]+30.0000)/5.00000)))
    algebraic[50] = (algebraic[49]*(power(constants[13]/(constants[13]+0.590000), 3.00000))*(states[0]+81.9000))/(1.00000+exp((1.39300*(states[0]+81.9000+3.60000))/constants[53]))
    algebraic[58] = exp((-constants[21]*states[0])/(2.00000*constants[53]))
    algebraic[53] = 1.00000+(constants[32]/constants[31])*(1.00000+exp((constants[22]*states[0])/constants[53]))+constants[33]/constants[29]+(power(constants[33], 2.00000))/(constants[29]*constants[30])+(power(constants[33], 3.00000))/(constants[29]*constants[30]*constants[28])
    algebraic[55] = (((power(constants[33], 2.00000))/(constants[29]*constants[30])+(power(constants[33], 3.00000))/(constants[29]*constants[30]*constants[28]))*exp((-constants[21]*states[0])/(2.00000*constants[53])))/algebraic[53]
    algebraic[56] = ((constants[32]/constants[31])*exp((-constants[22]*states[0])/constants[53]))/algebraic[53]
    algebraic[54] = exp((constants[21]*states[0])/(2.00000*constants[53]))
    algebraic[61] = algebraic[58]*constants[54]*(algebraic[55]+algebraic[56])+algebraic[56]*algebraic[54]*(constants[58]+algebraic[58])
    algebraic[57] = 1.00000+(states[5]/constants[23])*(1.00000+exp((-constants[20]*states[0])/constants[53])+constants[18]/constants[27])+constants[18]/constants[24]+(power(constants[18], 2.00000))/(constants[24]*constants[25])+(power(constants[18], 3.00000))/(constants[24]*constants[25]*constants[26])
    algebraic[60] = ((states[5]/constants[23])*exp((-constants[20]*states[0])/constants[53]))/algebraic[57]
    algebraic[59] = (((power(constants[18], 2.00000))/(constants[24]*constants[25])+(power(constants[18], 3.00000))/(constants[24]*constants[25]*constants[26]))*exp((constants[21]*states[0])/(2.00000*constants[53])))/algebraic[57]
    algebraic[62] = algebraic[54]*constants[58]*(algebraic[59]+algebraic[60])+algebraic[58]*algebraic[60]*(constants[54]+algebraic[54])
    algebraic[63] = algebraic[59]*constants[58]*(algebraic[55]+algebraic[56])+algebraic[60]*algebraic[55]*(constants[58]+algebraic[58])
    algebraic[64] = algebraic[55]*constants[54]*(algebraic[59]+algebraic[60])+algebraic[59]*algebraic[56]*(constants[54]+algebraic[54])
    algebraic[65] = (constants[19]*(algebraic[62]*algebraic[56]-algebraic[61]*algebraic[60]))/(algebraic[61]+algebraic[62]+algebraic[63]+algebraic[64])
    algebraic[52] = (constants[17]*(power(constants[18]/(5.64000+constants[18]), 3.00000))*(power(constants[13]/(0.621000+constants[13]), 2.00000))*1.60000)/(1.50000+exp(-(states[0]+60.0000)/40.0000))
    algebraic[51] = constants[15]*(states[0]-constants[16])
    algebraic[17] = custom_piecewise([greater_equal(voi , constants[4]) & less_equal(voi , constants[5]) & less_equal((voi-constants[4])-floor((voi-constants[4])/constants[6])*constants[6] , constants[7]), constants[8] , True, 0.00000])
    rates[0] = -(algebraic[66]+algebraic[71]+algebraic[67]+algebraic[42]+algebraic[34]+algebraic[68]+algebraic[50]+algebraic[65]+algebraic[52]+algebraic[51]+algebraic[70]+algebraic[17])/constants[3]
    algebraic[73] = 5.00000/(1.00000+constants[51]/states[19])
    algebraic[74] = (states[20]-states[21])/constants[52]
    rates[20] = algebraic[73]-(algebraic[74]*constants[60])/constants[59]
    algebraic[75] = 88800.0*states[19]*(1.00000-states[22])-446.000*states[22]
    rates[22] = algebraic[75]
    algebraic[77] = 227700.*states[19]*((1.00000-states[23])-states[24])-7.51000*states[23]
    rates[23] = algebraic[77]
    algebraic[76] = (constants[50]*(states[21]-states[5]))/(1.00000+power(0.00120000/states[5], 2.00000))
    algebraic[80] = 534.000*states[21]*(1.00000-states[27])-445.000*states[27]
    rates[21] = (algebraic[74]-algebraic[76])-10.0000*algebraic[80]
    algebraic[72] = (states[5]-states[19])/4.00000e-05
    algebraic[78] = 227700.*states[19]*(1.00000-states[25])-542.000*states[25]
    rates[19] = (algebraic[72]*constants[61]-algebraic[73]*constants[59])/constants[62]-(0.0450000*algebraic[78]+0.0310000*algebraic[75]+0.0620000*algebraic[77])
    rates[25] = algebraic[78]
    algebraic[79] = 227700.*states[5]*(1.00000-states[26])-542.000*states[26]
    rates[26] = algebraic[79]
    rates[27] = algebraic[80]
    algebraic[81] = 0.00100000*(115.000*states[5]*(1.00000-states[28])-1000.00*states[28])
    rates[5] = (((-(algebraic[71]-2.00000*algebraic[65])/(2.00000*constants[2])+algebraic[76]*constants[60])/constants[61]-algebraic[72])-0.0450000*algebraic[79])-(0.0310000/1.20000)*algebraic[81]
    rates[28] = algebraic[81]
    return(rates)

def computeAlgebraic(constants, states, voi):
    algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
    states = array(states)
    voi = array(voi)
    algebraic[15] = 120.000/(1.00000+exp(-(states[0]+50.0000)/15.0000))
    algebraic[16] = 5.82000/(1.00000+exp(-(states[0]+50.0000)/15.0000))
    algebraic[18] = 2277.00*2.50000*((1.00000-states[23])-states[24])-751.000*states[24]
    algebraic[0] = 1.00000/(1.00000+exp(((states[0]+83.1900)-(-7.20000*(power(constants[10], 0.690000)))/(power(1.26000e-05, 0.690000)+power(constants[10], 0.690000)))/13.5600))
    algebraic[19] = 0.250000+2.00000*exp(-(power(states[0]+70.0000, 2.00000))/500.000)
    algebraic[1] = 1.00000/(1.00000+exp((states[0]+10.2200)/-8.50000))
    algebraic[20] = 1.00000/(17.0000*exp(0.0398000*states[0])+0.211000*exp(-0.0510000*states[0]))
    algebraic[2] = 1.00000/(1.00000+exp((states[0]+10.2200)/-8.50000))
    algebraic[21] = 0.335810+0.906730*exp(-(power(states[0]+10.0000, 2.00000))/988.050)
    algebraic[8] = 1.00000/(1.00000+exp((states[0]-(-24.0000+constants[40]))/6.31000))
    algebraic[27] = 0.0100000+0.153900*exp(-(power(states[0]+40.0000, 2.00000))/185.670)
    algebraic[9] = 1.00000/(1.00000+exp((states[0]-(-24.0000+constants[41]))/6.31000))
    algebraic[28] = 0.0600000+0.480760*2.25000*exp(-(power(states[0]--40.0000, 2.00000))/138.040)
    algebraic[29] = 0.000596000+0.00311800/(1.03700*exp(0.0900000*(states[0]+30.6100))+0.396000*exp(-0.120000*(states[0]+23.8400)))
    algebraic[10] = 1.00000/(1.00000+exp((states[0]-7.44000)/-16.4000))
    algebraic[30] = 0.0126600+4.72716/(1.00000+exp((states[0]+154.500)/23.9600))
    algebraic[11] = 1.00000/(1.00000+exp((states[0]+33.8000)/6.12000))
    algebraic[31] = 0.100000+4.00000*exp(-(power(states[0]+65.0000, 2.00000))/500.000)
    algebraic[12] = 1.00000/(1.00000+exp((states[0]+33.8000)/6.12000))
    algebraic[4] = states[0]+44.4000
    algebraic[23] = custom_piecewise([less(fabs(algebraic[4]) , constants[35]), (-460.000*-12.6730)/exp(algebraic[4]/-12.6730) , True, (-460.000*algebraic[4])/(exp(algebraic[4]/-12.6730)-1.00000)])
    algebraic[36] = 18400.0*exp(algebraic[4]/-12.6730)
    algebraic[3] = (1.00000/(1.00000+exp((states[0]+4.90000)/15.1400)))*(1.00000-0.300000*exp(-(power(states[0], 2.00000))/500.000))
    algebraic[22] = 92.0100*exp(-0.0183000*states[0])
    algebraic[35] = 603.600*exp(0.00942000*states[0])
    algebraic[43] = 1.00000/(algebraic[22]+algebraic[35])
    algebraic[5] = 44.9000*exp((states[0]+66.9000)/-5.57000)
    algebraic[24] = 1491.00/(1.00000+323.300*exp((states[0]+94.6000)/-12.9000))
    algebraic[37] = algebraic[5]/(algebraic[5]+algebraic[24])
    algebraic[44] = 0.0300000/(1.00000+exp((states[0]+40.0000)/6.00000))+0.000350000
    algebraic[6] = 44.9000*exp((states[0]+66.9000)/-5.57000)
    algebraic[25] = 1491.00/(1.00000+323.300*exp((states[0]+94.6000)/-12.9000))
    algebraic[38] = algebraic[6]/(algebraic[6]+algebraic[25])
    algebraic[45] = 0.120000/(1.00000+exp((states[0]+60.0000)/2.00000))+0.00295000
    algebraic[39] = 1.00000/(1.00000+exp((states[0]-(-3.20000+constants[38]))/constants[39]))
    algebraic[7] = (-26.1200*(states[0]+35.0000))/(exp((states[0]+35.0000)/-2.50000)-1.00000)+(-78.1100*states[0])/(exp(-0.208000*states[0])-1.00000)
    algebraic[26] = (10.5200*(states[0]-5.00000))/(exp(0.400000*(states[0]-5.00000))-1.00000)
    algebraic[46] = 1.00000/(algebraic[7]+algebraic[26])
    algebraic[32] = 1.00000/(0.150000*exp(-states[0]/11.0000)+0.200000*exp(-states[0]/700.000))
    algebraic[40] = 1.00000/(16.0000*exp(states[0]/8.00000)+15.0000*exp(states[0]/50.0000))
    algebraic[47] = 0.00100000/(algebraic[32]+algebraic[40])
    algebraic[13] = 1.00000/(1.00000+exp((states[0]--49.1000)/-8.98000))
    algebraic[14] = 0.150400/(3100.00*exp(states[0]/13.0000)+700.000*exp(states[0]/70.0000))
    algebraic[33] = 0.150400/(95.0000*exp(-states[0]/10.0000)+50.0000*exp(-states[0]/700.000))+0.000229000/(1.00000+exp(-states[0]/5.00000))
    algebraic[48] = 0.00100000/(algebraic[14]+algebraic[33])
    algebraic[41] = algebraic[14]/(algebraic[14]+algebraic[33])
    algebraic[66] = (((constants[34]*(power(states[6], 3.00000))*(0.635000*states[7]+0.365000*states[8])*constants[33]*states[0]*constants[2])/constants[53])*(exp((states[0]-constants[56])/constants[53])-1.00000))/(exp(states[0]/constants[53])-1.00000)
    algebraic[69] = (constants[45]*states[17]*states[18]*(power(constants[10], 1.50000)))/(power(constants[46], 1.50000)+power(constants[10], 1.50000))
    algebraic[70] = (((algebraic[69]*constants[13])/(10.0000+constants[13]))*(states[0]-constants[55]))/(1.00000+exp(((states[0]-constants[55])-140.000)/(2.50000*constants[53])))
    algebraic[71] = constants[36]*states[9]*(0.675000*states[10]+0.325000*states[11])*(states[0]-constants[37])*(1.00000-((algebraic[70]*constants[10])/(9.00000e-05+constants[10]))/1.00000)
    algebraic[67] = constants[42]*states[12]*(0.450000*states[13]+0.550000*states[14])*(states[0]-constants[57])
    algebraic[42] = constants[11]*(0.900000*states[2]+0.100000*states[3])*states[4]*(states[0]-constants[55])
    algebraic[34] = states[1]*constants[9]*(states[0]--30.0000)
    algebraic[68] = constants[44]*states[15]*states[16]*(states[0]-constants[43])
    algebraic[49] = constants[14]*(0.500000+0.500000/(1.00000+exp((states[0]+30.0000)/5.00000)))
    algebraic[50] = (algebraic[49]*(power(constants[13]/(constants[13]+0.590000), 3.00000))*(states[0]+81.9000))/(1.00000+exp((1.39300*(states[0]+81.9000+3.60000))/constants[53]))
    algebraic[58] = exp((-constants[21]*states[0])/(2.00000*constants[53]))
    algebraic[53] = 1.00000+(constants[32]/constants[31])*(1.00000+exp((constants[22]*states[0])/constants[53]))+constants[33]/constants[29]+(power(constants[33], 2.00000))/(constants[29]*constants[30])+(power(constants[33], 3.00000))/(constants[29]*constants[30]*constants[28])
    algebraic[55] = (((power(constants[33], 2.00000))/(constants[29]*constants[30])+(power(constants[33], 3.00000))/(constants[29]*constants[30]*constants[28]))*exp((-constants[21]*states[0])/(2.00000*constants[53])))/algebraic[53]
    algebraic[56] = ((constants[32]/constants[31])*exp((-constants[22]*states[0])/constants[53]))/algebraic[53]
    algebraic[54] = exp((constants[21]*states[0])/(2.00000*constants[53]))
    algebraic[61] = algebraic[58]*constants[54]*(algebraic[55]+algebraic[56])+algebraic[56]*algebraic[54]*(constants[58]+algebraic[58])
    algebraic[57] = 1.00000+(states[5]/constants[23])*(1.00000+exp((-constants[20]*states[0])/constants[53])+constants[18]/constants[27])+constants[18]/constants[24]+(power(constants[18], 2.00000))/(constants[24]*constants[25])+(power(constants[18], 3.00000))/(constants[24]*constants[25]*constants[26])
    algebraic[60] = ((states[5]/constants[23])*exp((-constants[20]*states[0])/constants[53]))/algebraic[57]
    algebraic[59] = (((power(constants[18], 2.00000))/(constants[24]*constants[25])+(power(constants[18], 3.00000))/(constants[24]*constants[25]*constants[26]))*exp((constants[21]*states[0])/(2.00000*constants[53])))/algebraic[57]
    algebraic[62] = algebraic[54]*constants[58]*(algebraic[59]+algebraic[60])+algebraic[58]*algebraic[60]*(constants[54]+algebraic[54])
    algebraic[63] = algebraic[59]*constants[58]*(algebraic[55]+algebraic[56])+algebraic[60]*algebraic[55]*(constants[58]+algebraic[58])
    algebraic[64] = algebraic[55]*constants[54]*(algebraic[59]+algebraic[60])+algebraic[59]*algebraic[56]*(constants[54]+algebraic[54])
    algebraic[65] = (constants[19]*(algebraic[62]*algebraic[56]-algebraic[61]*algebraic[60]))/(algebraic[61]+algebraic[62]+algebraic[63]+algebraic[64])
    algebraic[52] = (constants[17]*(power(constants[18]/(5.64000+constants[18]), 3.00000))*(power(constants[13]/(0.621000+constants[13]), 2.00000))*1.60000)/(1.50000+exp(-(states[0]+60.0000)/40.0000))
    algebraic[51] = constants[15]*(states[0]-constants[16])
    algebraic[17] = custom_piecewise([greater_equal(voi , constants[4]) & less_equal(voi , constants[5]) & less_equal((voi-constants[4])-floor((voi-constants[4])/constants[6])*constants[6] , constants[7]), constants[8] , True, 0.00000])
    algebraic[73] = 5.00000/(1.00000+constants[51]/states[19])
    algebraic[74] = (states[20]-states[21])/constants[52]
    algebraic[75] = 88800.0*states[19]*(1.00000-states[22])-446.000*states[22]
    algebraic[77] = 227700.*states[19]*((1.00000-states[23])-states[24])-7.51000*states[23]
    algebraic[76] = (constants[50]*(states[21]-states[5]))/(1.00000+power(0.00120000/states[5], 2.00000))
    algebraic[80] = 534.000*states[21]*(1.00000-states[27])-445.000*states[27]
    algebraic[72] = (states[5]-states[19])/4.00000e-05
    algebraic[78] = 227700.*states[19]*(1.00000-states[25])-542.000*states[25]
    algebraic[79] = 227700.*states[5]*(1.00000-states[26])-542.000*states[26]
    algebraic[81] = 0.00100000*(115.000*states[5]*(1.00000-states[28])-1000.00*states[28])
    return algebraic

def custom_piecewise(cases):
    """Compute result of a piecewise function"""
    return select(cases[0::2],cases[1::2])

def solve_model():
    """Solve model with ODE solver"""
    from scipy.integrate import ode
    # Initialise constants and state variables
    (init_states, constants) = initConsts()

    # Set timespan to solve over
    voi = linspace(0, 10, 500)

    # Construct ODE object to solve
    r = ode(computeRates)
    r.set_integrator('vode', method='bdf', atol=1e-06, rtol=1e-06, max_step=1)
    r.set_initial_value(init_states, voi[0])
    r.set_f_params(constants)

    # Solve model
    states = array([[0.0] * len(voi)] * sizeStates)
    states[:,0] = init_states
    for (i,t) in enumerate(voi[1:]):
        if r.successful():
            r.integrate(t)
            states[:,i+1] = r.y
        else:
            break

    # Compute algebraic variables
    algebraic = computeAlgebraic(constants, states, voi)
    return (voi, states, algebraic)

def plot_model(voi, states, algebraic):
    """Plot variables against variable of integration"""
    import pylab
    (legend_states, legend_algebraic, legend_voi, legend_constants) = createLegends()
    pylab.figure(1)
    pylab.plot(voi,vstack((states,algebraic)).T)
    pylab.xlabel(legend_voi)
    pylab.legend(legend_states + legend_algebraic, loc='best')
    pylab.show()

if __name__ == "__main__":
    (voi, states, algebraic) = solve_model()
    plot_model(voi, states, algebraic)