function [VOI, STATES, ALGEBRAIC, CONSTANTS] = mainFunction()
    % This is the "main function".  In Matlab, things work best if you rename this function to match the filename.
   [VOI, STATES, ALGEBRAIC, CONSTANTS] = solveModel();
end

function [algebraicVariableCount] = getAlgebraicVariableCount() 
    % Used later when setting a global variable with the number of algebraic variables.
    % Note: This is not the "main method".  
    algebraicVariableCount =119;
end
% There are a total of 31 entries in each of the rate and state variable arrays.
% There are a total of 130 entries in the constant variable array.
%

function [VOI, STATES, ALGEBRAIC, CONSTANTS] = solveModel()
    % Create ALGEBRAIC of correct size
    global algebraicVariableCount;  algebraicVariableCount = getAlgebraicVariableCount();
    % Initialise constants and state variables
    [INIT_STATES, CONSTANTS] = initConsts;

    % Set timespan to solve over 
    tspan = [0, 10];

    % Set numerical accuracy options for ODE solver
    options = odeset('RelTol', 1e-06, 'AbsTol', 1e-06, 'MaxStep', 1);

    % Solve model with ODE solver
    [VOI, STATES] = ode15s(@(VOI, STATES)computeRates(VOI, STATES, CONSTANTS), tspan, INIT_STATES, options);

    % Compute algebraic variables
    [RATES, ALGEBRAIC] = computeRates(VOI, STATES, CONSTANTS);
    ALGEBRAIC = computeAlgebraic(ALGEBRAIC, CONSTANTS, STATES, VOI);

    % Plot state variables against variable of integration
    [LEGEND_STATES, LEGEND_ALGEBRAIC, LEGEND_VOI, LEGEND_CONSTANTS] = createLegends();
    figure();
    plot(VOI, STATES);
    xlabel(LEGEND_VOI);
    l = legend(LEGEND_STATES);
    set(l,'Interpreter','none');
end

function [LEGEND_STATES, LEGEND_ALGEBRAIC, LEGEND_VOI, LEGEND_CONSTANTS] = createLegends()
    LEGEND_STATES = ''; LEGEND_ALGEBRAIC = ''; LEGEND_VOI = ''; LEGEND_CONSTANTS = '';
    LEGEND_VOI = strpad('time in component environment (second)');
    LEGEND_CONSTANTS(:,1) = strpad('L_totmax in component b1_AR_Gs_parameters (uM)');
    LEGEND_CONSTANTS(:,2) = strpad('sum_b1_AR in component b1_AR_Gs_parameters (uM)');
    LEGEND_CONSTANTS(:,3) = strpad('Gs_tot in component b1_AR_Gs_parameters (uM)');
    LEGEND_CONSTANTS(:,4) = strpad('Kl in component b1_AR_Gs_parameters (uM)');
    LEGEND_CONSTANTS(:,5) = strpad('Kr in component b1_AR_Gs_parameters (uM)');
    LEGEND_CONSTANTS(:,6) = strpad('Kc in component b1_AR_Gs_parameters (uM)');
    LEGEND_CONSTANTS(:,7) = strpad('k_bar_kp in component b1_AR_Gs_parameters (per_sec)');
    LEGEND_CONSTANTS(:,8) = strpad('k_bar_km in component b1_AR_Gs_parameters (per_sec)');
    LEGEND_CONSTANTS(:,9) = strpad('k_p_kap in component b1_AR_Gs_parameters (per_uM_per_sec)');
    LEGEND_CONSTANTS(:,10) = strpad('k_p_kam in component b1_AR_Gs_parameters (per_sec)');
    LEGEND_CONSTANTS(:,11) = strpad('k_g_act in component b1_AR_Gs_parameters (per_sec)');
    LEGEND_CONSTANTS(:,12) = strpad('k_hyd in component b1_AR_Gs_parameters (per_sec)');
    LEGEND_CONSTANTS(:,13) = strpad('k_reassoc in component b1_AR_Gs_parameters (per_uM_per_sec)');
    LEGEND_CONSTANTS(:,14) = strpad('AC_tot in component cAMP_parameters (uM)');
    LEGEND_CONSTANTS(:,15) = strpad('ATP in component cAMP_parameters (uM)');
    LEGEND_CONSTANTS(:,16) = strpad('PDE_tot in component cAMP_parameters (uM)');
    LEGEND_CONSTANTS(:,17) = strpad('IBMX_tot in component cAMP_parameters (uM)');
    LEGEND_CONSTANTS(:,18) = strpad('Fsk_tot in component cAMP_parameters (uM)');
    LEGEND_CONSTANTS(:,19) = strpad('k_ac_basal in component cAMP_parameters (per_sec)');
    LEGEND_CONSTANTS(:,20) = strpad('k_ac_gsa in component cAMP_parameters (per_sec)');
    LEGEND_CONSTANTS(:,21) = strpad('k_ac_fsk in component cAMP_parameters (per_sec)');
    LEGEND_CONSTANTS(:,22) = strpad('k_pde in component cAMP_parameters (per_sec)');
    LEGEND_CONSTANTS(:,23) = strpad('Km_basal in component cAMP_parameters (uM)');
    LEGEND_CONSTANTS(:,24) = strpad('Km_gsa in component cAMP_parameters (uM)');
    LEGEND_CONSTANTS(:,25) = strpad('Km_fsk in component cAMP_parameters (uM)');
    LEGEND_CONSTANTS(:,26) = strpad('Km_pde in component cAMP_parameters (uM)');
    LEGEND_CONSTANTS(:,27) = strpad('K_gsa in component cAMP_parameters (uM)');
    LEGEND_CONSTANTS(:,28) = strpad('K_fsk in component cAMP_parameters (uM)');
    LEGEND_CONSTANTS(:,29) = strpad('Ki_ibmx in component cAMP_parameters (uM)');
    LEGEND_CONSTANTS(:,30) = strpad('PKAI_tot in component PKA_parameters (uM)');
    LEGEND_CONSTANTS(:,31) = strpad('PKAII_tot in component PKA_parameters (uM)');
    LEGEND_CONSTANTS(:,32) = strpad('PKI_tot in component PKA_parameters (uM)');
    LEGEND_CONSTANTS(:,33) = strpad('K_a in component PKA_parameters (uM)');
    LEGEND_CONSTANTS(:,34) = strpad('K_b in component PKA_parameters (uM)');
    LEGEND_CONSTANTS(:,35) = strpad('K_d in component PKA_parameters (uM)');
    LEGEND_CONSTANTS(:,36) = strpad('Ki_pki in component PKA_parameters (uM)');
    LEGEND_CONSTANTS(:,37) = strpad('PLB_tot in component PLB_parameters (uM)');
    LEGEND_CONSTANTS(:,38) = strpad('PP1_tot in component PLB_parameters (uM)');
    LEGEND_CONSTANTS(:,39) = strpad('Inhib1_tot in component PLB_parameters (uM)');
    LEGEND_CONSTANTS(:,40) = strpad('k_pka_plb in component PLB_parameters (per_sec)');
    LEGEND_CONSTANTS(:,41) = strpad('Km_pka_plb in component PLB_parameters (uM)');
    LEGEND_CONSTANTS(:,42) = strpad('k_pp1_plb in component PLB_parameters (per_sec)');
    LEGEND_CONSTANTS(:,43) = strpad('Km_pp1_plb in component PLB_parameters (uM)');
    LEGEND_CONSTANTS(:,44) = strpad('k_pka_i1 in component PLB_parameters (per_sec)');
    LEGEND_CONSTANTS(:,45) = strpad('Km_pka_i1 in component PLB_parameters (uM)');
    LEGEND_CONSTANTS(:,46) = strpad('Vmax_pp2a_i1 in component PLB_parameters (uM_per_sec)');
    LEGEND_CONSTANTS(:,47) = strpad('Km_pp2a_i1 in component PLB_parameters (uM)');
    LEGEND_CONSTANTS(:,48) = strpad('Ki_inhib1 in component PLB_parameters (uM)');
    LEGEND_CONSTANTS(:,49) = strpad('epsilon in component LCC_parameters (dimensionless)');
    LEGEND_CONSTANTS(:,50) = strpad('LCC_tot in component LCC_parameters (uM)');
    LEGEND_CONSTANTS(:,51) = strpad('PP1_lcc_tot in component LCC_parameters (uM)');
    LEGEND_CONSTANTS(:,52) = strpad('PP2A_lcc_tot in component LCC_parameters (uM)');
    LEGEND_CONSTANTS(:,53) = strpad('k_pka_lcc in component LCC_parameters (per_sec)');
    LEGEND_CONSTANTS(:,54) = strpad('Km_pka_lcc in component LCC_parameters (uM)');
    LEGEND_CONSTANTS(:,55) = strpad('k_pp1_lcc in component LCC_parameters (per_sec)');
    LEGEND_CONSTANTS(:,56) = strpad('Km_pp1_lcc in component LCC_parameters (uM)');
    LEGEND_CONSTANTS(:,57) = strpad('k_pp2a_lcc in component LCC_parameters (per_sec)');
    LEGEND_CONSTANTS(:,58) = strpad('Km_pp2a_lcc in component LCC_parameters (uM)');
    LEGEND_CONSTANTS(:,59) = strpad('V_myo in component EC_Coupling_Parameters (uL)');
    LEGEND_CONSTANTS(:,60) = strpad('V_nsr in component EC_Coupling_Parameters (uL)');
    LEGEND_CONSTANTS(:,61) = strpad('V_jsr in component EC_Coupling_Parameters (uL)');
    LEGEND_CONSTANTS(:,62) = strpad('A_Cap in component EC_Coupling_Parameters (cm2)');
    LEGEND_CONSTANTS(:,63) = strpad('Temp in component EC_Coupling_Parameters (kelvin)');
    LEGEND_CONSTANTS(:,64) = strpad('Na_ext in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,65) = strpad('K_ext in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,66) = strpad('Ca_ext in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,67) = strpad('G_Na in component EC_Coupling_Parameters (mS_per_uF)');
    LEGEND_CONSTANTS(:,68) = strpad('G_to in component EC_Coupling_Parameters (mS_per_uF)');
    LEGEND_CONSTANTS(:,69) = strpad('G_ss in component EC_Coupling_Parameters (mS_per_uF)');
    LEGEND_CONSTANTS(:,70) = strpad('G_Ki_bar in component EC_Coupling_Parameters (mS_per_uF)');
    LEGEND_CONSTANTS(:,71) = strpad('G_Kp in component EC_Coupling_Parameters (mS_per_uF)');
    LEGEND_CONSTANTS(:,72) = strpad('f in component EC_Coupling_Parameters (per_sec)');
    LEGEND_CONSTANTS(:,73) = strpad('g in component EC_Coupling_Parameters (per_sec)');
    LEGEND_CONSTANTS(:,74) = strpad('gamma_o in component EC_Coupling_Parameters (per_mM_per_sec)');
    LEGEND_CONSTANTS(:,75) = strpad('omega in component EC_Coupling_Parameters (per_sec)');
    LEGEND_CONSTANTS(:,76) = strpad('p_Ca in component EC_Coupling_Parameters (cm_per_sec)');
    LEGEND_CONSTANTS(:,77) = strpad('p_K in component EC_Coupling_Parameters (cm_per_sec)');
    LEGEND_CONSTANTS(:,78) = strpad('N_lcc in component EC_Coupling_Parameters (dimensionless)');
    LEGEND_CONSTANTS(:,79) = strpad('I_Ca05 in component EC_Coupling_Parameters (uA_per_uF)');
    LEGEND_CONSTANTS(:,80) = strpad('k_NaCa in component EC_Coupling_Parameters (uA_per_uF)');
    LEGEND_CONSTANTS(:,81) = strpad('Km_Na in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,82) = strpad('Km_Ca in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,83) = strpad('k_sat in component EC_Coupling_Parameters (dimensionless)');
    LEGEND_CONSTANTS(:,84) = strpad('eta in component EC_Coupling_Parameters (dimensionless)');
    LEGEND_CONSTANTS(:,85) = strpad('I_bar_NaK in component EC_Coupling_Parameters (uA_per_uF)');
    LEGEND_CONSTANTS(:,86) = strpad('Km_Nai in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,87) = strpad('Km_Ko in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,88) = strpad('I_bar_PCa in component EC_Coupling_Parameters (uA_per_uF)');
    LEGEND_CONSTANTS(:,89) = strpad('Km_PCa in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,90) = strpad('G_CaB in component EC_Coupling_Parameters (uA_per_uF)');
    LEGEND_CONSTANTS(:,91) = strpad('G_NaB in component EC_Coupling_Parameters (uA_per_uF)');
    LEGEND_CONSTANTS(:,92) = strpad('Pns in component EC_Coupling_Parameters (dimensionless)');
    LEGEND_CONSTANTS(:,93) = strpad('Km_NS in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,94) = strpad('I_up_bar in component EC_Coupling_Parameters (mM_per_sec)');
    LEGEND_CONSTANTS(:,95) = strpad('Km_up0 in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,96) = strpad('NSR_bar in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,97) = strpad('tau_on in component EC_Coupling_Parameters (second)');
    LEGEND_CONSTANTS(:,98) = strpad('tau_off in component EC_Coupling_Parameters (second)');
    LEGEND_CONSTANTS(:,99) = strpad('G_max_rel in component EC_Coupling_Parameters (mM_per_sec)');
    LEGEND_CONSTANTS(:,100) = strpad('d_Cai_th in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,101) = strpad('Km_rel in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,102) = strpad('CSQN_th in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,103) = strpad('CSQN_bar in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,104) = strpad('Km_CSQN in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,105) = strpad('tau_tr in component EC_Coupling_Parameters (second)');
    LEGEND_CONSTANTS(:,106) = strpad('TRPN_bar in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,107) = strpad('CMDN_bar in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,108) = strpad('INDO_bar in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,109) = strpad('Km_TRPN in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,110) = strpad('Km_CMDN in component EC_Coupling_Parameters (mM)');
    LEGEND_CONSTANTS(:,111) = strpad('Km_INDO in component EC_Coupling_Parameters (mM)');
    LEGEND_ALGEBRAIC(:,63) = strpad('LR in component b1_AR_module (uM)');
    LEGEND_ALGEBRAIC(:,64) = strpad('LRG in component b1_AR_module (uM)');
    LEGEND_ALGEBRAIC(:,65) = strpad('RG in component b1_AR_module (uM)');
    LEGEND_ALGEBRAIC(:,81) = strpad('BARK_DESENS in component b1_AR_module (uM_per_sec)');
    LEGEND_ALGEBRAIC(:,12) = strpad('BARK_RESENS in component b1_AR_module (uM_per_sec)');
    LEGEND_ALGEBRAIC(:,89) = strpad('PKA_DESENS in component b1_AR_module (uM_per_sec)');
    LEGEND_ALGEBRAIC(:,24) = strpad('PKA_RESENS in component b1_AR_module (uM_per_sec)');
    LEGEND_ALGEBRAIC(:,82) = strpad('G_ACT in component b1_AR_module (uM_per_sec)');
    LEGEND_ALGEBRAIC(:,26) = strpad('HYD in component b1_AR_module (uM_per_sec)');
    LEGEND_ALGEBRAIC(:,28) = strpad('REASSOC in component b1_AR_module (per_sec)');
    LEGEND_ALGEBRAIC(:,66) = strpad('L in component b1_AR_module (uM)');
    LEGEND_ALGEBRAIC(:,67) = strpad('R in component b1_AR_module (uM)');
    LEGEND_ALGEBRAIC(:,68) = strpad('Gs in component b1_AR_module (uM)');
    LEGEND_STATES(:,1) = strpad('b1_AR_tot in component b1_AR_module (uM)');
    LEGEND_STATES(:,2) = strpad('b1_AR_d in component b1_AR_module (uM)');
    LEGEND_STATES(:,3) = strpad('b1_AR_p in component b1_AR_module (uM)');
    LEGEND_STATES(:,4) = strpad('Gs_agtp_tot in component b1_AR_module (uM)');
    LEGEND_STATES(:,5) = strpad('Gs_agdp in component b1_AR_module (uM)');
    LEGEND_STATES(:,6) = strpad('Gs_bg in component b1_AR_module (uM)');
    LEGEND_ALGEBRAIC(:,69) = strpad('PKAC_I in component PKA_module (uM)');
    LEGEND_ALGEBRAIC(:,70) = strpad('cAMP in component PKA_module (uM)');
    LEGEND_ALGEBRAIC(:,43) = strpad('Gsa_GTP in component cAMP_module (uM)');
    LEGEND_ALGEBRAIC(:,51) = strpad('Fsk in component cAMP_module (uM)');
    LEGEND_ALGEBRAIC(:,44) = strpad('AC in component cAMP_module (uM)');
    LEGEND_CONSTANTS(:,128) = strpad('PDE in component cAMP_module (uM)');
    LEGEND_CONSTANTS(:,129) = strpad('IBMX in component cAMP_module (uM)');
    LEGEND_STATES(:,7) = strpad('cAMP_tot in component cAMP_module (uM)');
    LEGEND_ALGEBRAIC(:,45) = strpad('Gsa_GTP_AC in component cAMP_module (uM)');
    LEGEND_ALGEBRAIC(:,52) = strpad('Fsk_AC in component cAMP_module (uM)');
    LEGEND_ALGEBRAIC(:,49) = strpad('AC_ACT_GSA in component cAMP_module (uM)');
    LEGEND_ALGEBRAIC(:,47) = strpad('AC_ACT_BASAL in component cAMP_module (uM)');
    LEGEND_ALGEBRAIC(:,54) = strpad('AC_ACT_FSK in component cAMP_module (uM)');
    LEGEND_ALGEBRAIC(:,83) = strpad('PDE_ACT in component cAMP_module (uM)');
    LEGEND_CONSTANTS(:,130) = strpad('PDE_IBMX in component cAMP_module (uM)');
    LEGEND_ALGEBRAIC(:,71) = strpad('PKI in component PKA_module (uM)');
    LEGEND_ALGEBRAIC(:,72) = strpad('A2RC_I in component PKA_module (uM)');
    LEGEND_ALGEBRAIC(:,73) = strpad('A2R_I in component PKA_module (uM)');
    LEGEND_ALGEBRAIC(:,74) = strpad('A2RC_II in component PKA_module (uM)');
    LEGEND_ALGEBRAIC(:,75) = strpad('A2R_II in component PKA_module (uM)');
    LEGEND_ALGEBRAIC(:,76) = strpad('ARC_I in component PKA_module (uM)');
    LEGEND_ALGEBRAIC(:,77) = strpad('ARC_II in component PKA_module (uM)');
    LEGEND_ALGEBRAIC(:,78) = strpad('PKA_temp in component PKA_module (uM)');
    LEGEND_ALGEBRAIC(:,79) = strpad('PKAC_II in component PKA_module (uM)');
    LEGEND_ALGEBRAIC(:,30) = strpad('PLB in component PLB_module (uM)');
    LEGEND_ALGEBRAIC(:,84) = strpad('PLB_PHOSPH in component PLB_module (uM_per_sec)');
    LEGEND_ALGEBRAIC(:,61) = strpad('PLB_DEPHOSPH in component PLB_module (uM_per_sec)');
    LEGEND_ALGEBRAIC(:,31) = strpad('Inhib1 in component PLB_module (uM)');
    LEGEND_ALGEBRAIC(:,56) = strpad('Inhib1p_PP1 in component PLB_module (uM)');
    LEGEND_ALGEBRAIC(:,85) = strpad('Inhib1_PHOSPH in component PLB_module (uM_per_sec)');
    LEGEND_ALGEBRAIC(:,33) = strpad('Inhib1_DEPHOSPH in component PLB_module (uM_per_sec)');
    LEGEND_STATES(:,8) = strpad('PLB_p in component PLB_module (uM)');
    LEGEND_STATES(:,9) = strpad('Inhib1p_tot in component PLB_module (uM)');
    LEGEND_ALGEBRAIC(:,57) = strpad('Inhib1p in component PLB_module (uM)');
    LEGEND_ALGEBRAIC(:,58) = strpad('PP1 in component PLB_module (uM)');
    LEGEND_ALGEBRAIC(:,1) = strpad('frac_PLB_p in component PLB_module (dimensionless)');
    LEGEND_ALGEBRAIC(:,32) = strpad('frac_PLB in component PLB_module (dimensionless)');
    LEGEND_CONSTANTS(:,120) = strpad('frac_PLB_o in component PLB_module (dimensionless)');
    LEGEND_ALGEBRAIC(:,35) = strpad('LCCa in component LCC_module (uM)');
    LEGEND_ALGEBRAIC(:,86) = strpad('LCCa_PHOSPH in component LCC_module (uM_per_sec)');
    LEGEND_ALGEBRAIC(:,37) = strpad('LCCa_DEPHOSPH in component LCC_module (uM_per_sec)');
    LEGEND_ALGEBRAIC(:,39) = strpad('LCCb in component LCC_module (uM)');
    LEGEND_ALGEBRAIC(:,87) = strpad('LCCb_PHOSPH in component LCC_module (uM_per_sec)');
    LEGEND_ALGEBRAIC(:,41) = strpad('LCCb_DEPHOSPH in component LCC_module (uM_per_sec)');
    LEGEND_STATES(:,10) = strpad('LCCa_p in component LCC_module (uM)');
    LEGEND_STATES(:,11) = strpad('LCCb_p in component LCC_module (uM)');
    LEGEND_ALGEBRAIC(:,2) = strpad('frac_LCCa_p in component LCC_module (dimensionless)');
    LEGEND_CONSTANTS(:,121) = strpad('frac_LCCa_po in component LCC_module (dimensionless)');
    LEGEND_ALGEBRAIC(:,34) = strpad('frac_LCCb_p in component LCC_module (dimensionless)');
    LEGEND_CONSTANTS(:,122) = strpad('frac_LCCb_po in component LCC_module (dimensionless)');
    LEGEND_ALGEBRAIC(:,36) = strpad('E_Na in component Nernst_Potentials (mV)');
    LEGEND_ALGEBRAIC(:,38) = strpad('E_K in component Nernst_Potentials (mV)');
    LEGEND_ALGEBRAIC(:,40) = strpad('E_Ca in component Nernst_Potentials (mV)');
    LEGEND_CONSTANTS(:,123) = strpad('E_Cl in component Nernst_Potentials (mV)');
    LEGEND_CONSTANTS(:,112) = strpad('R in component Nernst_Potentials (joules_per_mole_kelvin)');
    LEGEND_CONSTANTS(:,113) = strpad('Frdy in component Nernst_Potentials (coulombs_per_mole)');
    LEGEND_CONSTANTS(:,124) = strpad('FoRT in component Nernst_Potentials (per_mV)');
    LEGEND_CONSTANTS(:,114) = strpad('z_Na in component Nernst_Potentials (dimensionless)');
    LEGEND_CONSTANTS(:,115) = strpad('z_K in component Nernst_Potentials (dimensionless)');
    LEGEND_CONSTANTS(:,116) = strpad('z_Ca in component Nernst_Potentials (dimensionless)');
    LEGEND_STATES(:,12) = strpad('Na_i in component Ion_Concentrations_and_Membrane_Potential (mM)');
    LEGEND_STATES(:,13) = strpad('K_i in component Ion_Concentrations_and_Membrane_Potential (mM)');
    LEGEND_STATES(:,14) = strpad('Ca_i in component Ion_Concentrations_and_Membrane_Potential (mM)');
    LEGEND_ALGEBRAIC(:,3) = strpad('am in component Fast_Na_Current (per_sec)');
    LEGEND_ALGEBRAIC(:,14) = strpad('bm in component Fast_Na_Current (per_sec)');
    LEGEND_ALGEBRAIC(:,4) = strpad('ah in component Fast_Na_Current (per_sec)');
    LEGEND_ALGEBRAIC(:,5) = strpad('aj in component Fast_Na_Current (per_sec)');
    LEGEND_ALGEBRAIC(:,15) = strpad('bh in component Fast_Na_Current (per_sec)');
    LEGEND_ALGEBRAIC(:,16) = strpad('bj in component Fast_Na_Current (per_sec)');
    LEGEND_STATES(:,15) = strpad('m in component Fast_Na_Current (dimensionless)');
    LEGEND_STATES(:,16) = strpad('h in component Fast_Na_Current (dimensionless)');
    LEGEND_STATES(:,17) = strpad('j in component Fast_Na_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,42) = strpad('I_Na in component Fast_Na_Current (uA_per_uF)');
    LEGEND_STATES(:,18) = strpad('V_m in component Ion_Concentrations_and_Membrane_Potential (mV)');
    LEGEND_ALGEBRAIC(:,6) = strpad('a_lcc in component L_Type_Calcium_Current (per_sec)');
    LEGEND_ALGEBRAIC(:,17) = strpad('b_lcc in component L_Type_Calcium_Current (per_sec)');
    LEGEND_ALGEBRAIC(:,19) = strpad('f_lcc in component L_Type_Calcium_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,7) = strpad('y_lcc_inf in component L_Type_Calcium_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,18) = strpad('tau_y_lcc in component L_Type_Calcium_Current (second)');
    LEGEND_ALGEBRAIC(:,25) = strpad('gamma in component L_Type_Calcium_Current (per_mM_per_sec)');
    LEGEND_ALGEBRAIC(:,27) = strpad('v_gamma in component L_Type_Calcium_Current (per_mM_per_sec)');
    LEGEND_ALGEBRAIC(:,29) = strpad('v_omega in component L_Type_Calcium_Current (dimensionless)');
    LEGEND_STATES(:,19) = strpad('v in component L_Type_Calcium_Current (dimensionless)');
    LEGEND_STATES(:,20) = strpad('w in component L_Type_Calcium_Current (dimensionless)');
    LEGEND_STATES(:,21) = strpad('x in component L_Type_Calcium_Current (dimensionless)');
    LEGEND_STATES(:,22) = strpad('y in component L_Type_Calcium_Current (dimensionless)');
    LEGEND_STATES(:,23) = strpad('z in component L_Type_Calcium_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,46) = strpad('i_bar_Ca in component L_Type_Calcium_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,48) = strpad('i_bar_K in component L_Type_Calcium_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,50) = strpad('f_avail in component L_Type_Calcium_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,53) = strpad('I_Ca in component L_Type_Calcium_Current (uA_per_uF)');
    LEGEND_ALGEBRAIC(:,55) = strpad('I_CaK in component L_Type_Calcium_Current (uA_per_uF)');
    LEGEND_ALGEBRAIC(:,59) = strpad('I_Ca_tot in component L_Type_Calcium_Current (uA_per_uF)');
    LEGEND_ALGEBRAIC(:,8) = strpad('r_toss in component Transient_Outward_K_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,9) = strpad('s_toss in component Transient_Outward_K_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,20) = strpad('tau_r_to in component Transient_Outward_K_Current (second)');
    LEGEND_ALGEBRAIC(:,21) = strpad('tau_s_to in component Transient_Outward_K_Current (second)');
    LEGEND_ALGEBRAIC(:,22) = strpad('tau_ss_to in component Transient_Outward_K_Current (second)');
    LEGEND_STATES(:,24) = strpad('r_to in component Transient_Outward_K_Current (dimensionless)');
    LEGEND_STATES(:,25) = strpad('s_to in component Transient_Outward_K_Current (dimensionless)');
    LEGEND_STATES(:,26) = strpad('ss_to in component Transient_Outward_K_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,60) = strpad('I_to in component Transient_Outward_K_Current (uA_per_uF)');
    LEGEND_ALGEBRAIC(:,10) = strpad('r_ss_inf in component Steady_State_K_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,23) = strpad('tau_r_ss in component Steady_State_K_Current (second)');
    LEGEND_ALGEBRAIC(:,11) = strpad('s_ss_inf in component Steady_State_K_Current (dimensionless)');
    LEGEND_CONSTANTS(:,125) = strpad('tau_s_ss in component Steady_State_K_Current (second)');
    LEGEND_STATES(:,27) = strpad('r_ss in component Steady_State_K_Current (dimensionless)');
    LEGEND_STATES(:,28) = strpad('s_ss in component Steady_State_K_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,62) = strpad('I_ss in component Steady_State_K_Current (uA_per_uF)');
    LEGEND_ALGEBRAIC(:,80) = strpad('a_Ki in component Time_Independent_K_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,88) = strpad('b_Ki in component Time_Independent_K_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,90) = strpad('Ki_ss in component Time_Independent_K_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,91) = strpad('I_Ki in component Time_Independent_K_Current (uA_per_uF)');
    LEGEND_ALGEBRAIC(:,92) = strpad('Kp in component Plateau_K_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,93) = strpad('I_Kp in component Plateau_K_Current (uA_per_uF)');
    LEGEND_ALGEBRAIC(:,94) = strpad('s4 in component Na_Ca_Exchanger_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,95) = strpad('s5 in component Na_Ca_Exchanger_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,96) = strpad('I_NCX in component Na_Ca_Exchanger_Current (uA_per_uF)');
    LEGEND_CONSTANTS(:,126) = strpad('sigma in component Na_K_Pump_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,97) = strpad('f_NaK in component Na_K_Pump_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,98) = strpad('I_NaK in component Na_K_Pump_Current (uA_per_uF)');
    LEGEND_ALGEBRAIC(:,99) = strpad('I_PCa in component Sarcolemmal_Ca_Pump_Current (uA_per_uF)');
    LEGEND_ALGEBRAIC(:,100) = strpad('I_CaB in component Ca_Background_Current (uA_per_uF)');
    LEGEND_ALGEBRAIC(:,101) = strpad('I_NaB in component Na_Background_Current (uA_per_uF)');
    LEGEND_ALGEBRAIC(:,102) = strpad('I_Na_tot in component Total_Membrane_Currents (uA_per_uF)');
    LEGEND_ALGEBRAIC(:,103) = strpad('I_K_tot in component Total_Membrane_Currents (uA_per_uF)');
    LEGEND_ALGEBRAIC(:,104) = strpad('I_Ca_tot in component Total_Membrane_Currents (uA_per_uF)');
    LEGEND_ALGEBRAIC(:,105) = strpad('t_rel in component Calcium_Induced_Calcium_Release (second)');
    LEGEND_ALGEBRAIC(:,106) = strpad('ryr_on in component Calcium_Induced_Calcium_Release (dimensionless)');
    LEGEND_ALGEBRAIC(:,107) = strpad('ryr_off in component Calcium_Induced_Calcium_Release (dimensionless)');
    LEGEND_ALGEBRAIC(:,108) = strpad('g_rel in component Calcium_Induced_Calcium_Release (per_sec)');
    LEGEND_ALGEBRAIC(:,109) = strpad('I_rel in component Calcium_Induced_Calcium_Release (mM_per_sec)');
    LEGEND_STATES(:,29) = strpad('Ca_jsr in component Other_SR_Fluxes_and_Concentrations (mM)');
    LEGEND_STATES(:,30) = strpad('trel in component Ion_Concentrations_and_Membrane_Potential (second)');
    LEGEND_ALGEBRAIC(:,110) = strpad('Km_up in component Other_SR_Fluxes_and_Concentrations (mM)');
    LEGEND_ALGEBRAIC(:,111) = strpad('I_up in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)');
    LEGEND_ALGEBRAIC(:,112) = strpad('I_leak in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)');
    LEGEND_ALGEBRAIC(:,113) = strpad('I_tr in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)');
    LEGEND_ALGEBRAIC(:,115) = strpad('B_jsr in component Other_SR_Fluxes_and_Concentrations (dimensionless)');
    LEGEND_STATES(:,31) = strpad('Ca_nsr in component Other_SR_Fluxes_and_Concentrations (mM)');
    LEGEND_ALGEBRAIC(:,117) = strpad('SR_content in component Other_SR_Fluxes_and_Concentrations (uM)');
    LEGEND_ALGEBRAIC(:,114) = strpad('b_trpn in component Cytoplasmic_Calcium_Buffering (dimensionless)');
    LEGEND_ALGEBRAIC(:,116) = strpad('b_cmdn in component Cytoplasmic_Calcium_Buffering (dimensionless)');
    LEGEND_ALGEBRAIC(:,118) = strpad('b_indo in component Cytoplasmic_Calcium_Buffering (dimensionless)');
    LEGEND_ALGEBRAIC(:,119) = strpad('B_myo in component Cytoplasmic_Calcium_Buffering (dimensionless)');
    LEGEND_CONSTANTS(:,127) = strpad('I_app in component Ion_Concentrations_and_Membrane_Potential (uA_per_uF)');
    LEGEND_CONSTANTS(:,117) = strpad('V_hold in component Ion_Concentrations_and_Membrane_Potential (mV)');
    LEGEND_CONSTANTS(:,118) = strpad('V_test in component Ion_Concentrations_and_Membrane_Potential (mV)');
    LEGEND_ALGEBRAIC(:,13) = strpad('V_clamp in component Ion_Concentrations_and_Membrane_Potential (mV)');
    LEGEND_CONSTANTS(:,119) = strpad('R_clamp in component Ion_Concentrations_and_Membrane_Potential (dimensionless)');
    LEGEND_RATES(:,1) = strpad('d/dt b1_AR_tot in component b1_AR_module (uM)');
    LEGEND_RATES(:,2) = strpad('d/dt b1_AR_d in component b1_AR_module (uM)');
    LEGEND_RATES(:,3) = strpad('d/dt b1_AR_p in component b1_AR_module (uM)');
    LEGEND_RATES(:,4) = strpad('d/dt Gs_agtp_tot in component b1_AR_module (uM)');
    LEGEND_RATES(:,5) = strpad('d/dt Gs_agdp in component b1_AR_module (uM)');
    LEGEND_RATES(:,6) = strpad('d/dt Gs_bg in component b1_AR_module (uM)');
    LEGEND_RATES(:,7) = strpad('d/dt cAMP_tot in component cAMP_module (uM)');
    LEGEND_RATES(:,8) = strpad('d/dt PLB_p in component PLB_module (uM)');
    LEGEND_RATES(:,9) = strpad('d/dt Inhib1p_tot in component PLB_module (uM)');
    LEGEND_RATES(:,10) = strpad('d/dt LCCa_p in component LCC_module (uM)');
    LEGEND_RATES(:,11) = strpad('d/dt LCCb_p in component LCC_module (uM)');
    LEGEND_RATES(:,15) = strpad('d/dt m in component Fast_Na_Current (dimensionless)');
    LEGEND_RATES(:,16) = strpad('d/dt h in component Fast_Na_Current (dimensionless)');
    LEGEND_RATES(:,17) = strpad('d/dt j in component Fast_Na_Current (dimensionless)');
    LEGEND_RATES(:,19) = strpad('d/dt v in component L_Type_Calcium_Current (dimensionless)');
    LEGEND_RATES(:,20) = strpad('d/dt w in component L_Type_Calcium_Current (dimensionless)');
    LEGEND_RATES(:,21) = strpad('d/dt x in component L_Type_Calcium_Current (dimensionless)');
    LEGEND_RATES(:,22) = strpad('d/dt y in component L_Type_Calcium_Current (dimensionless)');
    LEGEND_RATES(:,23) = strpad('d/dt z in component L_Type_Calcium_Current (dimensionless)');
    LEGEND_RATES(:,24) = strpad('d/dt r_to in component Transient_Outward_K_Current (dimensionless)');
    LEGEND_RATES(:,25) = strpad('d/dt s_to in component Transient_Outward_K_Current (dimensionless)');
    LEGEND_RATES(:,26) = strpad('d/dt ss_to in component Transient_Outward_K_Current (dimensionless)');
    LEGEND_RATES(:,27) = strpad('d/dt r_ss in component Steady_State_K_Current (dimensionless)');
    LEGEND_RATES(:,28) = strpad('d/dt s_ss in component Steady_State_K_Current (dimensionless)');
    LEGEND_RATES(:,31) = strpad('d/dt Ca_nsr in component Other_SR_Fluxes_and_Concentrations (mM)');
    LEGEND_RATES(:,29) = strpad('d/dt Ca_jsr in component Other_SR_Fluxes_and_Concentrations (mM)');
    LEGEND_RATES(:,12) = strpad('d/dt Na_i in component Ion_Concentrations_and_Membrane_Potential (mM)');
    LEGEND_RATES(:,13) = strpad('d/dt K_i in component Ion_Concentrations_and_Membrane_Potential (mM)');
    LEGEND_RATES(:,14) = strpad('d/dt Ca_i in component Ion_Concentrations_and_Membrane_Potential (mM)');
    LEGEND_RATES(:,18) = strpad('d/dt V_m in component Ion_Concentrations_and_Membrane_Potential (mV)');
    LEGEND_RATES(:,30) = strpad('d/dt trel in component Ion_Concentrations_and_Membrane_Potential (second)');
    LEGEND_STATES  = LEGEND_STATES';
    LEGEND_ALGEBRAIC = LEGEND_ALGEBRAIC';
    LEGEND_RATES = LEGEND_RATES';
    LEGEND_CONSTANTS = LEGEND_CONSTANTS';
end

function [STATES, CONSTANTS] = initConsts()
    VOI = 0; CONSTANTS = []; STATES = []; ALGEBRAIC = [];
    CONSTANTS(:,1) = 1;
    CONSTANTS(:,2) = 0.0132;
    CONSTANTS(:,3) = 3.83;
    CONSTANTS(:,4) = 0.285;
    CONSTANTS(:,5) = 0.062;
    CONSTANTS(:,6) = 33;
    CONSTANTS(:,7) = 1.1e-3;
    CONSTANTS(:,8) = 2.2e-3;
    CONSTANTS(:,9) = 3.6e-3;
    CONSTANTS(:,10) = 2.2e-3;
    CONSTANTS(:,11) = 16;
    CONSTANTS(:,12) = 0.8;
    CONSTANTS(:,13) = 1.21e3;
    CONSTANTS(:,14) = 49.7e-3;
    CONSTANTS(:,15) = 5e3;
    CONSTANTS(:,16) = 38.9e-3;
    CONSTANTS(:,17) = 0;
    CONSTANTS(:,18) = 0;
    CONSTANTS(:,19) = 0.2;
    CONSTANTS(:,20) = 8.5;
    CONSTANTS(:,21) = 7.3;
    CONSTANTS(:,22) = 5;
    CONSTANTS(:,23) = 1.03e3;
    CONSTANTS(:,24) = 315;
    CONSTANTS(:,25) = 860;
    CONSTANTS(:,26) = 1.3;
    CONSTANTS(:,27) = 0.4;
    CONSTANTS(:,28) = 44;
    CONSTANTS(:,29) = 30;
    CONSTANTS(:,30) = 0.59;
    CONSTANTS(:,31) = 0.025;
    CONSTANTS(:,32) = 0.18;
    CONSTANTS(:,33) = 9.14;
    CONSTANTS(:,34) = 1.64;
    CONSTANTS(:,35) = 4.375;
    CONSTANTS(:,36) = 0.2e-3;
    CONSTANTS(:,37) = 106;
    CONSTANTS(:,38) = 0.89;
    CONSTANTS(:,39) = 0.3;
    CONSTANTS(:,40) = 54;
    CONSTANTS(:,41) = 21;
    CONSTANTS(:,42) = 8.5;
    CONSTANTS(:,43) = 7;
    CONSTANTS(:,44) = 60;
    CONSTANTS(:,45) = 1;
    CONSTANTS(:,46) = 14;
    CONSTANTS(:,47) = 1;
    CONSTANTS(:,48) = 1e-3;
    CONSTANTS(:,49) = 10;
    CONSTANTS(:,50) = 25e-3;
    CONSTANTS(:,51) = 25e-3;
    CONSTANTS(:,52) = 25e-3;
    CONSTANTS(:,53) = 54;
    CONSTANTS(:,54) = 21;
    CONSTANTS(:,55) = 8.52;
    CONSTANTS(:,56) = 3;
    CONSTANTS(:,57) = 10.1;
    CONSTANTS(:,58) = 3;
    CONSTANTS(:,59) = 20.8e-6;
    CONSTANTS(:,60) = 9.88e-7;
    CONSTANTS(:,61) = 9.3e-8;
    CONSTANTS(:,62) = 1.534e-4;
    CONSTANTS(:,63) = 310;
    CONSTANTS(:,64) = 140;
    CONSTANTS(:,65) = 5.4;
    CONSTANTS(:,66) = 1.8;
    CONSTANTS(:,67) = 8;
    CONSTANTS(:,68) = 0.35;
    CONSTANTS(:,69) = 0.07;
    CONSTANTS(:,70) = 0.24;
    CONSTANTS(:,71) = 0.008;
    CONSTANTS(:,72) = 300;
    CONSTANTS(:,73) = 2e3;
    CONSTANTS(:,74) = 5187.5;
    CONSTANTS(:,75) = 10;
    CONSTANTS(:,76) = 1.7469e-8;
    CONSTANTS(:,77) = 3.234e-11;
    CONSTANTS(:,78) = 3e5;
    CONSTANTS(:,79) = -0.458;
    CONSTANTS(:,80) = 1483;
    CONSTANTS(:,81) = 87.5;
    CONSTANTS(:,82) = 1.38;
    CONSTANTS(:,83) = 0.1;
    CONSTANTS(:,84) = 0.35;
    CONSTANTS(:,85) = 1.1;
    CONSTANTS(:,86) = 10;
    CONSTANTS(:,87) = 1.5;
    CONSTANTS(:,88) = 1.15;
    CONSTANTS(:,89) = 0.5e-3;
    CONSTANTS(:,90) = 2.8e-3;
    CONSTANTS(:,91) = 1.18e-3;
    CONSTANTS(:,92) = 0;
    CONSTANTS(:,93) = 1.2e-3;
    CONSTANTS(:,94) = 4.7;
    CONSTANTS(:,95) = 3e-4;
    CONSTANTS(:,96) = 15;
    CONSTANTS(:,97) = 2e-3;
    CONSTANTS(:,98) = 2e-3;
    CONSTANTS(:,99) = 60e3;
    CONSTANTS(:,100) = 0.18e-3;
    CONSTANTS(:,101) = 0.8e-3;
    CONSTANTS(:,102) = 8.75;
    CONSTANTS(:,103) = 15;
    CONSTANTS(:,104) = 0.8;
    CONSTANTS(:,105) = 5.7e-4;
    CONSTANTS(:,106) = 0.07;
    CONSTANTS(:,107) = 0.05;
    CONSTANTS(:,108) = 0.07;
    CONSTANTS(:,109) = 0.5128e-3;
    CONSTANTS(:,110) = 2.38e-3;
    CONSTANTS(:,111) = 8.44e-4;
    STATES(:,1) = 0.01205;
    STATES(:,2) = 0;
    STATES(:,3) = 1.154e-3;
    STATES(:,4) = 0.02505;
    STATES(:,5) = 6.446e-4;
    STATES(:,6) = 0.02569;
    STATES(:,7) = 0.8453;
    STATES(:,8) = 4.105;
    STATES(:,9) = 0.0526;
    STATES(:,10) = 5.103e-3;
    STATES(:,11) = 5.841e-3;
    CONSTANTS(:,112) = 8314;
    CONSTANTS(:,113) = 96485;
    CONSTANTS(:,114) = 1;
    CONSTANTS(:,115) = 1;
    CONSTANTS(:,116) = 2;
    STATES(:,12) = 16;
    STATES(:,13) = 145;
    STATES(:,14) = 1.58e-4;
    STATES(:,15) = 1.4e-3;
    STATES(:,16) = 0.99;
    STATES(:,17) = 0.99;
    STATES(:,18) = -85.66;
    STATES(:,19) = 0;
    STATES(:,20) = 0;
    STATES(:,21) = 0.13;
    STATES(:,22) = 0.96;
    STATES(:,23) = 0.92;
    STATES(:,24) = 1.4e-3;
    STATES(:,25) = 1;
    STATES(:,26) = 0.613;
    STATES(:,27) = 198e-3;
    STATES(:,28) = 0.43;
    STATES(:,29) = 1.92;
    STATES(:,30) = 0.9;
    STATES(:,31) = 1.92;
    CONSTANTS(:,117) = -40;
    CONSTANTS(:,118) = -10;
    CONSTANTS(:,119) = 0.02;
    CONSTANTS(:,120) = 0.961300;
    CONSTANTS(:,121) = 0.204100;
    CONSTANTS(:,122) = 0.233600;
    CONSTANTS(:,123) =  - 40.0000;
    CONSTANTS(:,124) = (CONSTANTS(:,113)./CONSTANTS(:,112))./CONSTANTS(:,63);
    CONSTANTS(:,125) = 2.10000;
    CONSTANTS(:,126) = (exp(CONSTANTS(:,64)./67.3000) - 1.00000)./7.00000;
    CONSTANTS(:,127) = 0.00000;
    [CONSTANTS, STATES, ALGEBRAIC] = rootfind_0(VOI, CONSTANTS, STATES, ALGEBRAIC);
    if (isempty(STATES)), warning('Initial values for states not set');, end
end

function [RATES, ALGEBRAIC] = computeRates(VOI, STATES, CONSTANTS)
    global algebraicVariableCount;
    statesSize = size(STATES);
    statesColumnCount = statesSize(2);
    if ( statesColumnCount == 1)
        STATES = STATES';
        ALGEBRAIC = zeros(1, algebraicVariableCount);
        utilOnes = 1;
    else
        statesRowCount = statesSize(1);
        ALGEBRAIC = zeros(statesRowCount, algebraicVariableCount);
        RATES = zeros(statesRowCount, statesColumnCount);
        utilOnes = ones(statesRowCount, 1);
    end
    ALGEBRAIC(:,11) = 1.00000./(1.00000+exp((STATES(:,18)+87.5000)./10.3000));
    RATES(:,28) = (ALGEBRAIC(:,11) - STATES(:,28))./CONSTANTS(:,125);
    ALGEBRAIC(:,3) = ( 0.320000.*(STATES(:,18)+47.1300))./(1.00000 - exp(  - 0.100000.*(STATES(:,18)+47.1300)));
    ALGEBRAIC(:,14) =  0.0800000.*exp( - STATES(:,18)./11.0000);
    RATES(:,15) =  1000.00.*( ALGEBRAIC(:,3).*(1.00000 - STATES(:,15)) -  ALGEBRAIC(:,14).*STATES(:,15));
    ALGEBRAIC(:,4) = piecewise({STATES(:,18)>= - 40.0000, 0.00000 },  0.135000.*exp((80.0000+STATES(:,18))./ - 6.80000));
    ALGEBRAIC(:,15) = piecewise({STATES(:,18)>= - 40.0000, 1.00000./( 0.130000.*(1.00000+exp( - (STATES(:,18)+10.6600)./11.1000))) },  3.56000.*exp( 0.0790000.*STATES(:,18))+ 310000..*exp( 0.350000.*STATES(:,18)));
    RATES(:,16) =  1000.00.*( ALGEBRAIC(:,4).*(1.00000 - STATES(:,16)) -  ALGEBRAIC(:,15).*STATES(:,16));
    ALGEBRAIC(:,5) = piecewise({STATES(:,18)>= - 40.0000, 0.00000 }, ( (  - 127140..*exp( 0.244400.*STATES(:,18)) -  3.47400e-05.*exp(  - 0.0439100.*STATES(:,18))).*(STATES(:,18)+37.7800))./(1.00000+exp( 0.311000.*(STATES(:,18)+79.2300))));
    ALGEBRAIC(:,16) = piecewise({STATES(:,18)>= - 40.0000, ( 0.300000.*exp(  - 2.57500e-07.*STATES(:,18)))./(1.00000+exp(  - 0.100000.*(STATES(:,18)+32.0000))) }, ( 0.121200.*exp(  - 0.0105200.*STATES(:,18)))./(1.00000+exp(  - 0.137800.*(STATES(:,18)+40.1400))));
    RATES(:,17) =  1000.00.*( ALGEBRAIC(:,5).*(1.00000 - STATES(:,17)) -  ALGEBRAIC(:,16).*STATES(:,17));
    ALGEBRAIC(:,6) =  400.000.*exp((STATES(:,18)+2.00000)./10.0000);
    ALGEBRAIC(:,17) =  50.0000.*exp((  - 1.00000.*(STATES(:,18)+2.00000))./13.0000);
    RATES(:,19) =  ALGEBRAIC(:,6).*(1.00000 - STATES(:,19)) -  ALGEBRAIC(:,17).*STATES(:,19);
    RATES(:,20) =  2.00000.*ALGEBRAIC(:,6).*(1.00000 - STATES(:,20)) -  (ALGEBRAIC(:,17)./2.00000).*STATES(:,20);
    ALGEBRAIC(:,2) = STATES(:,10)./CONSTANTS(:,50);
    ALGEBRAIC(:,19) =  CONSTANTS(:,72).*(( 0.375000.*ALGEBRAIC(:,2))./CONSTANTS(:,121)+0.625000);
    RATES(:,21) =  ALGEBRAIC(:,19).*(1.00000 - STATES(:,21)) -  CONSTANTS(:,73).*STATES(:,21);
    ALGEBRAIC(:,7) = 1.00000./(1.00000+exp((STATES(:,18)+55.0000)./7.50000))+0.100000./(1.00000+exp(( - STATES(:,18)+21.0000)./6.00000));
    ALGEBRAIC(:,18) = 0.0200000+0.300000./(1.00000+exp((STATES(:,18)+30.0000)./9.50000));
    RATES(:,22) = (ALGEBRAIC(:,7) - STATES(:,22))./ALGEBRAIC(:,18);
    ALGEBRAIC(:,8) = 1.00000./(1.00000+exp((STATES(:,18)+10.6000)./ - 11.4200));
    ALGEBRAIC(:,20) = 1.00000./( 45.1600.*exp( 0.0357700.*(STATES(:,18)+50.0000))+ 98.9000.*exp(  - 0.100000.*(STATES(:,18)+38.0000)));
    RATES(:,24) = (ALGEBRAIC(:,8) - STATES(:,24))./ALGEBRAIC(:,20);
    ALGEBRAIC(:,9) = 1.00000./(1.00000+exp((STATES(:,18)+43.5000)./6.88410));
    ALGEBRAIC(:,21) =  0.350000.*exp(  - 1.00000.*power((STATES(:,18)+70.0000)./15.0000, 2.00000))+0.0350000;
    RATES(:,25) = (ALGEBRAIC(:,9) - STATES(:,25))./ALGEBRAIC(:,21);
    ALGEBRAIC(:,22) =  3.70000.*exp(  - 1.00000.*power((STATES(:,18)+70.0000)./30.0000, 2.00000))+0.0350000;
    RATES(:,26) = (ALGEBRAIC(:,9) - STATES(:,26))./ALGEBRAIC(:,22);
    ALGEBRAIC(:,10) = 1.00000./(1.00000+exp((STATES(:,18)+11.5000)./ - 11.8200));
    ALGEBRAIC(:,23) = 10.0000./( 45.1600.*exp( 0.0357700.*(STATES(:,18)+50.0000))+ 98.9000.*exp(  - 0.100000.*(STATES(:,18)+38.0000)));
    RATES(:,27) = (ALGEBRAIC(:,10) - STATES(:,27))./ALGEBRAIC(:,23);
    ALGEBRAIC(:,26) =  CONSTANTS(:,12).*STATES(:,4);
    ALGEBRAIC(:,28) =  CONSTANTS(:,13).*STATES(:,5).*STATES(:,6);
    RATES(:,5) = ALGEBRAIC(:,26) - ALGEBRAIC(:,28);
    ALGEBRAIC(:,25) =  CONSTANTS(:,74).*STATES(:,14);
    ALGEBRAIC(:,27) =  ALGEBRAIC(:,25).*(power(1.00000 - STATES(:,19), 4.00000)+ 2.00000.*STATES(:,19).*power(1.00000 - STATES(:,19), 3.00000)+ 4.00000.*power(STATES(:,19), 2.00000).*power(1.00000 - STATES(:,19), 2.00000)+ 8.00000.*power(STATES(:,19), 3.00000).*(1.00000 - STATES(:,19))+ 16.0000.*power(STATES(:,19), 4.00000).*(1.00000 - ALGEBRAIC(:,19)./CONSTANTS(:,73)));
    ALGEBRAIC(:,29) =  CONSTANTS(:,75).*(power(1.00000 - STATES(:,20), 4.00000)+ 0.500000.*STATES(:,20).*power(1.00000 - STATES(:,20), 3.00000)+ 0.250000.*power(STATES(:,20), 2.00000).*power(1.00000 - STATES(:,20), 2.00000)+ 0.125000.*power(STATES(:,20), 3.00000).*(1.00000 - STATES(:,20))+ 0.0625000.*power(STATES(:,20), 4.00000));
    RATES(:,23) =  ALGEBRAIC(:,29).*(1.00000 - STATES(:,23)) -  ALGEBRAIC(:,27).*STATES(:,23);
    [CONSTANTS, STATES, ALGEBRAIC] = rootfind_1(VOI, CONSTANTS, STATES, ALGEBRAIC);
    ALGEBRAIC(:,81) =  CONSTANTS(:,7).*(ALGEBRAIC(:,63)+ALGEBRAIC(:,64));
    ALGEBRAIC(:,12) =  CONSTANTS(:,8).*STATES(:,2);
    RATES(:,2) = ALGEBRAIC(:,81) - ALGEBRAIC(:,12);
    ALGEBRAIC(:,82) =  CONSTANTS(:,11).*(ALGEBRAIC(:,65)+ALGEBRAIC(:,64));
    RATES(:,4) = ALGEBRAIC(:,82) - ALGEBRAIC(:,26);
    RATES(:,6) = ALGEBRAIC(:,82) - ALGEBRAIC(:,28);
    [CONSTANTS, STATES, ALGEBRAIC] = rootfind_2(VOI, CONSTANTS, STATES, ALGEBRAIC);
    ALGEBRAIC(:,49) = ( CONSTANTS(:,20).*ALGEBRAIC(:,45).*CONSTANTS(:,15))./(CONSTANTS(:,24)+CONSTANTS(:,15));
    ALGEBRAIC(:,47) = ( CONSTANTS(:,19).*ALGEBRAIC(:,44).*CONSTANTS(:,15))./(CONSTANTS(:,23)+CONSTANTS(:,15));
    [CONSTANTS, STATES, ALGEBRAIC] = rootfind_3(VOI, CONSTANTS, STATES, ALGEBRAIC);
    ALGEBRAIC(:,54) = ( CONSTANTS(:,21).*ALGEBRAIC(:,52).*CONSTANTS(:,15))./(CONSTANTS(:,25)+CONSTANTS(:,15));
    ALGEBRAIC(:,83) = ( CONSTANTS(:,22).*CONSTANTS(:,128).*ALGEBRAIC(:,70))./(CONSTANTS(:,26)+ALGEBRAIC(:,70));
    RATES(:,7) = (ALGEBRAIC(:,47)+ALGEBRAIC(:,49)+ALGEBRAIC(:,54)) - ALGEBRAIC(:,83);
    ALGEBRAIC(:,30) = CONSTANTS(:,37) - STATES(:,8);
    ALGEBRAIC(:,84) = ( CONSTANTS(:,40).*ALGEBRAIC(:,69).*ALGEBRAIC(:,30))./(CONSTANTS(:,41)+ALGEBRAIC(:,30));
    [CONSTANTS, STATES, ALGEBRAIC] = rootfind_4(VOI, CONSTANTS, STATES, ALGEBRAIC);
    ALGEBRAIC(:,61) = ( CONSTANTS(:,42).*ALGEBRAIC(:,58).*STATES(:,8))./(CONSTANTS(:,43)+STATES(:,8));
    RATES(:,8) = ALGEBRAIC(:,84) - ALGEBRAIC(:,61);
    ALGEBRAIC(:,31) = CONSTANTS(:,39) - STATES(:,9);
    ALGEBRAIC(:,85) = ( CONSTANTS(:,44).*ALGEBRAIC(:,69).*ALGEBRAIC(:,31))./(CONSTANTS(:,45)+ALGEBRAIC(:,31));
    ALGEBRAIC(:,33) = ( CONSTANTS(:,46).*STATES(:,9))./(CONSTANTS(:,47)+STATES(:,9));
    RATES(:,9) = ALGEBRAIC(:,85) - ALGEBRAIC(:,33);
    ALGEBRAIC(:,35) = CONSTANTS(:,50) - STATES(:,10);
    ALGEBRAIC(:,86) = ( CONSTANTS(:,49).*CONSTANTS(:,53).*ALGEBRAIC(:,79).*ALGEBRAIC(:,35))./(CONSTANTS(:,54)+ CONSTANTS(:,49).*ALGEBRAIC(:,35));
    ALGEBRAIC(:,37) = ( CONSTANTS(:,49).*CONSTANTS(:,57).*CONSTANTS(:,52).*STATES(:,10))./(CONSTANTS(:,58)+ CONSTANTS(:,49).*STATES(:,10));
    RATES(:,10) = ALGEBRAIC(:,86) - ALGEBRAIC(:,37);
    ALGEBRAIC(:,39) = CONSTANTS(:,50) - STATES(:,11);
    ALGEBRAIC(:,87) = ( CONSTANTS(:,49).*CONSTANTS(:,53).*ALGEBRAIC(:,79).*ALGEBRAIC(:,39))./(CONSTANTS(:,54)+ CONSTANTS(:,49).*ALGEBRAIC(:,39));
    ALGEBRAIC(:,41) = ( CONSTANTS(:,49).*CONSTANTS(:,55).*CONSTANTS(:,51).*STATES(:,11))./(CONSTANTS(:,56)+ CONSTANTS(:,49).*STATES(:,11));
    RATES(:,11) = ALGEBRAIC(:,87) - ALGEBRAIC(:,41);
    ALGEBRAIC(:,89) =  CONSTANTS(:,9).*ALGEBRAIC(:,69).*STATES(:,1);
    ALGEBRAIC(:,24) =  CONSTANTS(:,10).*STATES(:,3);
    RATES(:,1) = (ALGEBRAIC(:,12) - ALGEBRAIC(:,81))+(ALGEBRAIC(:,24) - ALGEBRAIC(:,89));
    RATES(:,3) = ALGEBRAIC(:,89) - ALGEBRAIC(:,24);
    ALGEBRAIC(:,36) =  (1.00000./( CONSTANTS(:,124).*CONSTANTS(:,114))).*log(CONSTANTS(:,64)./STATES(:,12));
    ALGEBRAIC(:,42) =  CONSTANTS(:,67).*power(STATES(:,15), 3.00000).*STATES(:,16).*STATES(:,17).*(STATES(:,18) - ALGEBRAIC(:,36));
    ALGEBRAIC(:,94) =  exp( CONSTANTS(:,84).*STATES(:,18).*CONSTANTS(:,124)).*power(STATES(:,12), 3.00000).*CONSTANTS(:,66);
    ALGEBRAIC(:,95) =  exp( (CONSTANTS(:,84) - 1.00000).*STATES(:,18).*CONSTANTS(:,124)).*power(CONSTANTS(:,64), 3.00000).*STATES(:,14);
    ALGEBRAIC(:,96) =  (CONSTANTS(:,80)./( (power(CONSTANTS(:,81), 3.00000)+power(CONSTANTS(:,64), 3.00000)).*(CONSTANTS(:,82)+CONSTANTS(:,66)).*(1.00000+ CONSTANTS(:,83).*exp( (CONSTANTS(:,84) - 1.00000).*STATES(:,18).*CONSTANTS(:,124))))).*(ALGEBRAIC(:,94) - ALGEBRAIC(:,95));
    ALGEBRAIC(:,97) = 1.00000./(1.00000+ 0.124500.*exp(  - 0.100000.*STATES(:,18).*CONSTANTS(:,124))+ 0.0365000.*CONSTANTS(:,126).*exp(  - STATES(:,18).*CONSTANTS(:,124)));
    ALGEBRAIC(:,98) = ( (( CONSTANTS(:,85).*ALGEBRAIC(:,97))./(1.00000+power(CONSTANTS(:,86)./STATES(:,12), 1.50000))).*CONSTANTS(:,65))./(CONSTANTS(:,65)+CONSTANTS(:,87));
    ALGEBRAIC(:,101) =  CONSTANTS(:,91).*(STATES(:,18) - ALGEBRAIC(:,36));
    ALGEBRAIC(:,102) = ALGEBRAIC(:,42)+ALGEBRAIC(:,101)+ 3.00000.*ALGEBRAIC(:,96)+ 3.00000.*ALGEBRAIC(:,98);
    RATES(:,12) = (  - 1000.00.*ALGEBRAIC(:,102).*CONSTANTS(:,62))./( CONSTANTS(:,59).*CONSTANTS(:,114).*CONSTANTS(:,113));
    ALGEBRAIC(:,48) = ( CONSTANTS(:,77).*STATES(:,18).*CONSTANTS(:,113).*CONSTANTS(:,124).*( STATES(:,13).*exp( STATES(:,18).*CONSTANTS(:,124)) - CONSTANTS(:,65)))./(exp( STATES(:,18).*CONSTANTS(:,124)) - 1.00000);
    ALGEBRAIC(:,34) = STATES(:,11)./CONSTANTS(:,50);
    ALGEBRAIC(:,50) =  0.500000.*(( 0.400000.*ALGEBRAIC(:,34))./CONSTANTS(:,122)+0.600000);
    ALGEBRAIC(:,46) = ( CONSTANTS(:,76).*4.00000.*STATES(:,18).*CONSTANTS(:,113).*CONSTANTS(:,124).*( 0.00100000.*exp( 2.00000.*STATES(:,18).*CONSTANTS(:,124)) -  0.341000.*CONSTANTS(:,66)))./(exp( 2.00000.*STATES(:,18).*CONSTANTS(:,124)) - 1.00000);
    ALGEBRAIC(:,53) =  ALGEBRAIC(:,46).*CONSTANTS(:,78).*ALGEBRAIC(:,50).*power(STATES(:,19), 4.00000).*STATES(:,21).*STATES(:,22).*STATES(:,23);
    ALGEBRAIC(:,55) =  (ALGEBRAIC(:,48)./(1.00000+ALGEBRAIC(:,53)./CONSTANTS(:,79))).*CONSTANTS(:,78).*ALGEBRAIC(:,50).*power(STATES(:,19), 4.00000).*STATES(:,21).*STATES(:,22).*STATES(:,23);
    ALGEBRAIC(:,38) =  (1.00000./( CONSTANTS(:,124).*CONSTANTS(:,115))).*log(CONSTANTS(:,65)./STATES(:,13));
    ALGEBRAIC(:,60) =  CONSTANTS(:,68).*STATES(:,24).*( 0.886000.*STATES(:,25)+ 0.114000.*STATES(:,26)).*(STATES(:,18) - ALGEBRAIC(:,38));
    ALGEBRAIC(:,62) =  CONSTANTS(:,69).*STATES(:,27).*STATES(:,28).*(STATES(:,18) - ALGEBRAIC(:,38));
    ALGEBRAIC(:,80) = 1.02000./(1.00000+exp( 0.238500.*((STATES(:,18) - ALGEBRAIC(:,38)) - 59.2150)));
    ALGEBRAIC(:,88) = ( 0.491240.*exp( 0.0803200.*((STATES(:,18)+5.47600) - ALGEBRAIC(:,38)))+exp( 0.0617500.*((STATES(:,18) - ALGEBRAIC(:,38)) - 594.310)))./(1.00000+exp(  - 0.514300.*((STATES(:,18) - ALGEBRAIC(:,38))+4.75300)));
    ALGEBRAIC(:,90) = ALGEBRAIC(:,80)./(ALGEBRAIC(:,80)+ALGEBRAIC(:,88));
    ALGEBRAIC(:,91) =  CONSTANTS(:,70).*power((CONSTANTS(:,65)./5.40000), 1.0 ./ 2).*ALGEBRAIC(:,90).*(STATES(:,18) - ALGEBRAIC(:,38));
    ALGEBRAIC(:,92) = 1.00000./(1.00000+exp((7.48800 - STATES(:,18))./5.98000));
    ALGEBRAIC(:,93) =  CONSTANTS(:,71).*ALGEBRAIC(:,92).*(STATES(:,18) - ALGEBRAIC(:,38));
    ALGEBRAIC(:,103) = ((ALGEBRAIC(:,60)+ALGEBRAIC(:,62)+ALGEBRAIC(:,91)+ALGEBRAIC(:,93)) -  2.00000.*ALGEBRAIC(:,98))+ALGEBRAIC(:,55);
    RATES(:,13) = (  - 1000.00.*ALGEBRAIC(:,103).*CONSTANTS(:,62))./( CONSTANTS(:,59).*CONSTANTS(:,115).*CONSTANTS(:,113));
    ALGEBRAIC(:,99) = ( CONSTANTS(:,88).*STATES(:,14))./(CONSTANTS(:,89)+STATES(:,14));
    ALGEBRAIC(:,40) =  (1.00000./( CONSTANTS(:,124).*CONSTANTS(:,116))).*log(CONSTANTS(:,66)./STATES(:,14));
    ALGEBRAIC(:,100) =  CONSTANTS(:,90).*(STATES(:,18) - ALGEBRAIC(:,40));
    ALGEBRAIC(:,104) = (ALGEBRAIC(:,53)+ALGEBRAIC(:,100)+ALGEBRAIC(:,99)) -  2.00000.*ALGEBRAIC(:,96);
    RATES(:,18) =   - 1000.00.*((ALGEBRAIC(:,104)+ALGEBRAIC(:,103)+ALGEBRAIC(:,102)) - CONSTANTS(:,127));
    RATES(:,30) = piecewise({  - 1000.00.*((ALGEBRAIC(:,104)+ALGEBRAIC(:,103)+ALGEBRAIC(:,102)) - CONSTANTS(:,127))>30000.0, 1.00000 -  10000.0.*STATES(:,30) }, 1.00000);
    ALGEBRAIC(:,32) = ALGEBRAIC(:,30)./CONSTANTS(:,37);
    ALGEBRAIC(:,110) = ( CONSTANTS(:,95).*(1.00000+ 2.00000.*ALGEBRAIC(:,32)))./(1.00000+ 2.00000.*CONSTANTS(:,120));
    ALGEBRAIC(:,111) = ( CONSTANTS(:,94).*power(STATES(:,14), 2.00000))./(power(ALGEBRAIC(:,110), 2.00000)+power(STATES(:,14), 2.00000));
    ALGEBRAIC(:,112) = ( CONSTANTS(:,94).*STATES(:,31))./CONSTANTS(:,96);
    ALGEBRAIC(:,113) = (STATES(:,31) - STATES(:,29))./CONSTANTS(:,105);
    RATES(:,31) = (ALGEBRAIC(:,111) - ALGEBRAIC(:,112)) - ( ALGEBRAIC(:,113).*CONSTANTS(:,61))./CONSTANTS(:,60);
    ALGEBRAIC(:,105) = STATES(:,30)+0.00200000;
    ALGEBRAIC(:,106) = 1.00000 - exp( - ALGEBRAIC(:,105)./CONSTANTS(:,97));
    ALGEBRAIC(:,107) = exp( - ALGEBRAIC(:,105)./CONSTANTS(:,98));
    ALGEBRAIC(:,108) = CONSTANTS(:,99)./(1.00000+exp((ALGEBRAIC(:,104)+5.00000)./0.900000));
    ALGEBRAIC(:,109) =  ALGEBRAIC(:,108).*ALGEBRAIC(:,106).*ALGEBRAIC(:,107).*(STATES(:,29) - STATES(:,14));
    ALGEBRAIC(:,115) = 1.00000./(1.00000+( CONSTANTS(:,103).*CONSTANTS(:,104))./power(CONSTANTS(:,104)+STATES(:,29), 2.00000));
    RATES(:,29) =  ALGEBRAIC(:,115).*(ALGEBRAIC(:,113) - ALGEBRAIC(:,109));
    ALGEBRAIC(:,114) = ( CONSTANTS(:,106).*CONSTANTS(:,109))./power(CONSTANTS(:,109)+STATES(:,14), 2.00000);
    ALGEBRAIC(:,116) = ( CONSTANTS(:,107).*CONSTANTS(:,110))./power(CONSTANTS(:,110)+STATES(:,14), 2.00000);
    ALGEBRAIC(:,118) = ( CONSTANTS(:,108).*CONSTANTS(:,111))./power(CONSTANTS(:,111)+STATES(:,14), 2.00000);
    ALGEBRAIC(:,119) = 1.00000./(1.00000+ALGEBRAIC(:,116)+ALGEBRAIC(:,114)+ALGEBRAIC(:,114)+ALGEBRAIC(:,118));
    RATES(:,14) =   - ALGEBRAIC(:,119).*((( 1000.00.*ALGEBRAIC(:,104).*CONSTANTS(:,62))./( CONSTANTS(:,59).*CONSTANTS(:,116).*CONSTANTS(:,113))+( (ALGEBRAIC(:,111) - ALGEBRAIC(:,112)).*CONSTANTS(:,60))./CONSTANTS(:,59)) - ( ALGEBRAIC(:,109).*CONSTANTS(:,61))./CONSTANTS(:,59));
   RATES = RATES';
end

% Calculate algebraic variables
function ALGEBRAIC = computeAlgebraic(ALGEBRAIC, CONSTANTS, STATES, VOI)
    statesSize = size(STATES);
    statesColumnCount = statesSize(2);
    if ( statesColumnCount == 1)
        STATES = STATES';
        utilOnes = 1;
    else
        statesRowCount = statesSize(1);
        utilOnes = ones(statesRowCount, 1);
    end
    ALGEBRAIC(:,11) = 1.00000./(1.00000+exp((STATES(:,18)+87.5000)./10.3000));
    ALGEBRAIC(:,3) = ( 0.320000.*(STATES(:,18)+47.1300))./(1.00000 - exp(  - 0.100000.*(STATES(:,18)+47.1300)));
    ALGEBRAIC(:,14) =  0.0800000.*exp( - STATES(:,18)./11.0000);
    ALGEBRAIC(:,4) = piecewise({STATES(:,18)>= - 40.0000, 0.00000 },  0.135000.*exp((80.0000+STATES(:,18))./ - 6.80000));
    ALGEBRAIC(:,15) = piecewise({STATES(:,18)>= - 40.0000, 1.00000./( 0.130000.*(1.00000+exp( - (STATES(:,18)+10.6600)./11.1000))) },  3.56000.*exp( 0.0790000.*STATES(:,18))+ 310000..*exp( 0.350000.*STATES(:,18)));
    ALGEBRAIC(:,5) = piecewise({STATES(:,18)>= - 40.0000, 0.00000 }, ( (  - 127140..*exp( 0.244400.*STATES(:,18)) -  3.47400e-05.*exp(  - 0.0439100.*STATES(:,18))).*(STATES(:,18)+37.7800))./(1.00000+exp( 0.311000.*(STATES(:,18)+79.2300))));
    ALGEBRAIC(:,16) = piecewise({STATES(:,18)>= - 40.0000, ( 0.300000.*exp(  - 2.57500e-07.*STATES(:,18)))./(1.00000+exp(  - 0.100000.*(STATES(:,18)+32.0000))) }, ( 0.121200.*exp(  - 0.0105200.*STATES(:,18)))./(1.00000+exp(  - 0.137800.*(STATES(:,18)+40.1400))));
    ALGEBRAIC(:,6) =  400.000.*exp((STATES(:,18)+2.00000)./10.0000);
    ALGEBRAIC(:,17) =  50.0000.*exp((  - 1.00000.*(STATES(:,18)+2.00000))./13.0000);
    ALGEBRAIC(:,2) = STATES(:,10)./CONSTANTS(:,50);
    ALGEBRAIC(:,19) =  CONSTANTS(:,72).*(( 0.375000.*ALGEBRAIC(:,2))./CONSTANTS(:,121)+0.625000);
    ALGEBRAIC(:,7) = 1.00000./(1.00000+exp((STATES(:,18)+55.0000)./7.50000))+0.100000./(1.00000+exp(( - STATES(:,18)+21.0000)./6.00000));
    ALGEBRAIC(:,18) = 0.0200000+0.300000./(1.00000+exp((STATES(:,18)+30.0000)./9.50000));
    ALGEBRAIC(:,8) = 1.00000./(1.00000+exp((STATES(:,18)+10.6000)./ - 11.4200));
    ALGEBRAIC(:,20) = 1.00000./( 45.1600.*exp( 0.0357700.*(STATES(:,18)+50.0000))+ 98.9000.*exp(  - 0.100000.*(STATES(:,18)+38.0000)));
    ALGEBRAIC(:,9) = 1.00000./(1.00000+exp((STATES(:,18)+43.5000)./6.88410));
    ALGEBRAIC(:,21) =  0.350000.*exp(  - 1.00000.*power((STATES(:,18)+70.0000)./15.0000, 2.00000))+0.0350000;
    ALGEBRAIC(:,22) =  3.70000.*exp(  - 1.00000.*power((STATES(:,18)+70.0000)./30.0000, 2.00000))+0.0350000;
    ALGEBRAIC(:,10) = 1.00000./(1.00000+exp((STATES(:,18)+11.5000)./ - 11.8200));
    ALGEBRAIC(:,23) = 10.0000./( 45.1600.*exp( 0.0357700.*(STATES(:,18)+50.0000))+ 98.9000.*exp(  - 0.100000.*(STATES(:,18)+38.0000)));
    ALGEBRAIC(:,26) =  CONSTANTS(:,12).*STATES(:,4);
    ALGEBRAIC(:,28) =  CONSTANTS(:,13).*STATES(:,5).*STATES(:,6);
    ALGEBRAIC(:,25) =  CONSTANTS(:,74).*STATES(:,14);
    ALGEBRAIC(:,27) =  ALGEBRAIC(:,25).*(power(1.00000 - STATES(:,19), 4.00000)+ 2.00000.*STATES(:,19).*power(1.00000 - STATES(:,19), 3.00000)+ 4.00000.*power(STATES(:,19), 2.00000).*power(1.00000 - STATES(:,19), 2.00000)+ 8.00000.*power(STATES(:,19), 3.00000).*(1.00000 - STATES(:,19))+ 16.0000.*power(STATES(:,19), 4.00000).*(1.00000 - ALGEBRAIC(:,19)./CONSTANTS(:,73)));
    ALGEBRAIC(:,29) =  CONSTANTS(:,75).*(power(1.00000 - STATES(:,20), 4.00000)+ 0.500000.*STATES(:,20).*power(1.00000 - STATES(:,20), 3.00000)+ 0.250000.*power(STATES(:,20), 2.00000).*power(1.00000 - STATES(:,20), 2.00000)+ 0.125000.*power(STATES(:,20), 3.00000).*(1.00000 - STATES(:,20))+ 0.0625000.*power(STATES(:,20), 4.00000));
    ALGEBRAIC(:,81) =  CONSTANTS(:,7).*(ALGEBRAIC(:,63)+ALGEBRAIC(:,64));
    ALGEBRAIC(:,12) =  CONSTANTS(:,8).*STATES(:,2);
    ALGEBRAIC(:,82) =  CONSTANTS(:,11).*(ALGEBRAIC(:,65)+ALGEBRAIC(:,64));
    ALGEBRAIC(:,49) = ( CONSTANTS(:,20).*ALGEBRAIC(:,45).*CONSTANTS(:,15))./(CONSTANTS(:,24)+CONSTANTS(:,15));
    ALGEBRAIC(:,47) = ( CONSTANTS(:,19).*ALGEBRAIC(:,44).*CONSTANTS(:,15))./(CONSTANTS(:,23)+CONSTANTS(:,15));
    ALGEBRAIC(:,54) = ( CONSTANTS(:,21).*ALGEBRAIC(:,52).*CONSTANTS(:,15))./(CONSTANTS(:,25)+CONSTANTS(:,15));
    ALGEBRAIC(:,83) = ( CONSTANTS(:,22).*CONSTANTS(:,128).*ALGEBRAIC(:,70))./(CONSTANTS(:,26)+ALGEBRAIC(:,70));
    ALGEBRAIC(:,30) = CONSTANTS(:,37) - STATES(:,8);
    ALGEBRAIC(:,84) = ( CONSTANTS(:,40).*ALGEBRAIC(:,69).*ALGEBRAIC(:,30))./(CONSTANTS(:,41)+ALGEBRAIC(:,30));
    ALGEBRAIC(:,61) = ( CONSTANTS(:,42).*ALGEBRAIC(:,58).*STATES(:,8))./(CONSTANTS(:,43)+STATES(:,8));
    ALGEBRAIC(:,31) = CONSTANTS(:,39) - STATES(:,9);
    ALGEBRAIC(:,85) = ( CONSTANTS(:,44).*ALGEBRAIC(:,69).*ALGEBRAIC(:,31))./(CONSTANTS(:,45)+ALGEBRAIC(:,31));
    ALGEBRAIC(:,33) = ( CONSTANTS(:,46).*STATES(:,9))./(CONSTANTS(:,47)+STATES(:,9));
    ALGEBRAIC(:,35) = CONSTANTS(:,50) - STATES(:,10);
    ALGEBRAIC(:,86) = ( CONSTANTS(:,49).*CONSTANTS(:,53).*ALGEBRAIC(:,79).*ALGEBRAIC(:,35))./(CONSTANTS(:,54)+ CONSTANTS(:,49).*ALGEBRAIC(:,35));
    ALGEBRAIC(:,37) = ( CONSTANTS(:,49).*CONSTANTS(:,57).*CONSTANTS(:,52).*STATES(:,10))./(CONSTANTS(:,58)+ CONSTANTS(:,49).*STATES(:,10));
    ALGEBRAIC(:,39) = CONSTANTS(:,50) - STATES(:,11);
    ALGEBRAIC(:,87) = ( CONSTANTS(:,49).*CONSTANTS(:,53).*ALGEBRAIC(:,79).*ALGEBRAIC(:,39))./(CONSTANTS(:,54)+ CONSTANTS(:,49).*ALGEBRAIC(:,39));
    ALGEBRAIC(:,41) = ( CONSTANTS(:,49).*CONSTANTS(:,55).*CONSTANTS(:,51).*STATES(:,11))./(CONSTANTS(:,56)+ CONSTANTS(:,49).*STATES(:,11));
    ALGEBRAIC(:,89) =  CONSTANTS(:,9).*ALGEBRAIC(:,69).*STATES(:,1);
    ALGEBRAIC(:,24) =  CONSTANTS(:,10).*STATES(:,3);
    ALGEBRAIC(:,36) =  (1.00000./( CONSTANTS(:,124).*CONSTANTS(:,114))).*log(CONSTANTS(:,64)./STATES(:,12));
    ALGEBRAIC(:,42) =  CONSTANTS(:,67).*power(STATES(:,15), 3.00000).*STATES(:,16).*STATES(:,17).*(STATES(:,18) - ALGEBRAIC(:,36));
    ALGEBRAIC(:,94) =  exp( CONSTANTS(:,84).*STATES(:,18).*CONSTANTS(:,124)).*power(STATES(:,12), 3.00000).*CONSTANTS(:,66);
    ALGEBRAIC(:,95) =  exp( (CONSTANTS(:,84) - 1.00000).*STATES(:,18).*CONSTANTS(:,124)).*power(CONSTANTS(:,64), 3.00000).*STATES(:,14);
    ALGEBRAIC(:,96) =  (CONSTANTS(:,80)./( (power(CONSTANTS(:,81), 3.00000)+power(CONSTANTS(:,64), 3.00000)).*(CONSTANTS(:,82)+CONSTANTS(:,66)).*(1.00000+ CONSTANTS(:,83).*exp( (CONSTANTS(:,84) - 1.00000).*STATES(:,18).*CONSTANTS(:,124))))).*(ALGEBRAIC(:,94) - ALGEBRAIC(:,95));
    ALGEBRAIC(:,97) = 1.00000./(1.00000+ 0.124500.*exp(  - 0.100000.*STATES(:,18).*CONSTANTS(:,124))+ 0.0365000.*CONSTANTS(:,126).*exp(  - STATES(:,18).*CONSTANTS(:,124)));
    ALGEBRAIC(:,98) = ( (( CONSTANTS(:,85).*ALGEBRAIC(:,97))./(1.00000+power(CONSTANTS(:,86)./STATES(:,12), 1.50000))).*CONSTANTS(:,65))./(CONSTANTS(:,65)+CONSTANTS(:,87));
    ALGEBRAIC(:,101) =  CONSTANTS(:,91).*(STATES(:,18) - ALGEBRAIC(:,36));
    ALGEBRAIC(:,102) = ALGEBRAIC(:,42)+ALGEBRAIC(:,101)+ 3.00000.*ALGEBRAIC(:,96)+ 3.00000.*ALGEBRAIC(:,98);
    ALGEBRAIC(:,48) = ( CONSTANTS(:,77).*STATES(:,18).*CONSTANTS(:,113).*CONSTANTS(:,124).*( STATES(:,13).*exp( STATES(:,18).*CONSTANTS(:,124)) - CONSTANTS(:,65)))./(exp( STATES(:,18).*CONSTANTS(:,124)) - 1.00000);
    ALGEBRAIC(:,34) = STATES(:,11)./CONSTANTS(:,50);
    ALGEBRAIC(:,50) =  0.500000.*(( 0.400000.*ALGEBRAIC(:,34))./CONSTANTS(:,122)+0.600000);
    ALGEBRAIC(:,46) = ( CONSTANTS(:,76).*4.00000.*STATES(:,18).*CONSTANTS(:,113).*CONSTANTS(:,124).*( 0.00100000.*exp( 2.00000.*STATES(:,18).*CONSTANTS(:,124)) -  0.341000.*CONSTANTS(:,66)))./(exp( 2.00000.*STATES(:,18).*CONSTANTS(:,124)) - 1.00000);
    ALGEBRAIC(:,53) =  ALGEBRAIC(:,46).*CONSTANTS(:,78).*ALGEBRAIC(:,50).*power(STATES(:,19), 4.00000).*STATES(:,21).*STATES(:,22).*STATES(:,23);
    ALGEBRAIC(:,55) =  (ALGEBRAIC(:,48)./(1.00000+ALGEBRAIC(:,53)./CONSTANTS(:,79))).*CONSTANTS(:,78).*ALGEBRAIC(:,50).*power(STATES(:,19), 4.00000).*STATES(:,21).*STATES(:,22).*STATES(:,23);
    ALGEBRAIC(:,38) =  (1.00000./( CONSTANTS(:,124).*CONSTANTS(:,115))).*log(CONSTANTS(:,65)./STATES(:,13));
    ALGEBRAIC(:,60) =  CONSTANTS(:,68).*STATES(:,24).*( 0.886000.*STATES(:,25)+ 0.114000.*STATES(:,26)).*(STATES(:,18) - ALGEBRAIC(:,38));
    ALGEBRAIC(:,62) =  CONSTANTS(:,69).*STATES(:,27).*STATES(:,28).*(STATES(:,18) - ALGEBRAIC(:,38));
    ALGEBRAIC(:,80) = 1.02000./(1.00000+exp( 0.238500.*((STATES(:,18) - ALGEBRAIC(:,38)) - 59.2150)));
    ALGEBRAIC(:,88) = ( 0.491240.*exp( 0.0803200.*((STATES(:,18)+5.47600) - ALGEBRAIC(:,38)))+exp( 0.0617500.*((STATES(:,18) - ALGEBRAIC(:,38)) - 594.310)))./(1.00000+exp(  - 0.514300.*((STATES(:,18) - ALGEBRAIC(:,38))+4.75300)));
    ALGEBRAIC(:,90) = ALGEBRAIC(:,80)./(ALGEBRAIC(:,80)+ALGEBRAIC(:,88));
    ALGEBRAIC(:,91) =  CONSTANTS(:,70).*power((CONSTANTS(:,65)./5.40000), 1.0 ./ 2).*ALGEBRAIC(:,90).*(STATES(:,18) - ALGEBRAIC(:,38));
    ALGEBRAIC(:,92) = 1.00000./(1.00000+exp((7.48800 - STATES(:,18))./5.98000));
    ALGEBRAIC(:,93) =  CONSTANTS(:,71).*ALGEBRAIC(:,92).*(STATES(:,18) - ALGEBRAIC(:,38));
    ALGEBRAIC(:,103) = ((ALGEBRAIC(:,60)+ALGEBRAIC(:,62)+ALGEBRAIC(:,91)+ALGEBRAIC(:,93)) -  2.00000.*ALGEBRAIC(:,98))+ALGEBRAIC(:,55);
    ALGEBRAIC(:,99) = ( CONSTANTS(:,88).*STATES(:,14))./(CONSTANTS(:,89)+STATES(:,14));
    ALGEBRAIC(:,40) =  (1.00000./( CONSTANTS(:,124).*CONSTANTS(:,116))).*log(CONSTANTS(:,66)./STATES(:,14));
    ALGEBRAIC(:,100) =  CONSTANTS(:,90).*(STATES(:,18) - ALGEBRAIC(:,40));
    ALGEBRAIC(:,104) = (ALGEBRAIC(:,53)+ALGEBRAIC(:,100)+ALGEBRAIC(:,99)) -  2.00000.*ALGEBRAIC(:,96);
    ALGEBRAIC(:,32) = ALGEBRAIC(:,30)./CONSTANTS(:,37);
    ALGEBRAIC(:,110) = ( CONSTANTS(:,95).*(1.00000+ 2.00000.*ALGEBRAIC(:,32)))./(1.00000+ 2.00000.*CONSTANTS(:,120));
    ALGEBRAIC(:,111) = ( CONSTANTS(:,94).*power(STATES(:,14), 2.00000))./(power(ALGEBRAIC(:,110), 2.00000)+power(STATES(:,14), 2.00000));
    ALGEBRAIC(:,112) = ( CONSTANTS(:,94).*STATES(:,31))./CONSTANTS(:,96);
    ALGEBRAIC(:,113) = (STATES(:,31) - STATES(:,29))./CONSTANTS(:,105);
    ALGEBRAIC(:,105) = STATES(:,30)+0.00200000;
    ALGEBRAIC(:,106) = 1.00000 - exp( - ALGEBRAIC(:,105)./CONSTANTS(:,97));
    ALGEBRAIC(:,107) = exp( - ALGEBRAIC(:,105)./CONSTANTS(:,98));
    ALGEBRAIC(:,108) = CONSTANTS(:,99)./(1.00000+exp((ALGEBRAIC(:,104)+5.00000)./0.900000));
    ALGEBRAIC(:,109) =  ALGEBRAIC(:,108).*ALGEBRAIC(:,106).*ALGEBRAIC(:,107).*(STATES(:,29) - STATES(:,14));
    ALGEBRAIC(:,115) = 1.00000./(1.00000+( CONSTANTS(:,103).*CONSTANTS(:,104))./power(CONSTANTS(:,104)+STATES(:,29), 2.00000));
    ALGEBRAIC(:,114) = ( CONSTANTS(:,106).*CONSTANTS(:,109))./power(CONSTANTS(:,109)+STATES(:,14), 2.00000);
    ALGEBRAIC(:,116) = ( CONSTANTS(:,107).*CONSTANTS(:,110))./power(CONSTANTS(:,110)+STATES(:,14), 2.00000);
    ALGEBRAIC(:,118) = ( CONSTANTS(:,108).*CONSTANTS(:,111))./power(CONSTANTS(:,111)+STATES(:,14), 2.00000);
    ALGEBRAIC(:,119) = 1.00000./(1.00000+ALGEBRAIC(:,116)+ALGEBRAIC(:,114)+ALGEBRAIC(:,114)+ALGEBRAIC(:,118));
    ALGEBRAIC(:,1) = STATES(:,8)./CONSTANTS(:,37);
    ALGEBRAIC(:,13) = piecewise({VOI>59.1000&VOI<59.5000, CONSTANTS(:,118) }, CONSTANTS(:,117));
    ALGEBRAIC(:,59) = ALGEBRAIC(:,53)+ALGEBRAIC(:,55);
    ALGEBRAIC(:,117) =  1000.00.*(( (STATES(:,29)+STATES(:,29)./ALGEBRAIC(:,115)).*CONSTANTS(:,61))./CONSTANTS(:,59)+( STATES(:,31).*CONSTANTS(:,60))./CONSTANTS(:,59));
end

% Functions required for solving differential algebraic equation
function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_0(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN)
    ALGEBRAIC = ALGEBRAIC_IN;
    CONSTANTS = CONSTANTS_IN;
    STATES = STATES_IN;
    global initialGuess_0;
    if (length(initialGuess_0) ~= 3), initialGuess_0 = [0.0389,0,0.1];, end
    options = optimset('Display', 'off', 'TolX', 1E-6);
    if length(VOI) == 1
        residualfn = @(algebraicCandidate)residualSN_0(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES);
        soln = fsolve(residualfn, initialGuess_0, options);
        initialGuess_0 = soln;
        CONSTANTS(:,128) = soln(1);
        CONSTANTS(:,129) = soln(2);
        CONSTANTS(:,130) = soln(3);
    else
        SET_CONSTANTS(:,128) = logical(1);
        SET_CONSTANTS(:,129) = logical(1);
        SET_CONSTANTS(:,130) = logical(1);
        for i=1:length(VOI)
            residualfn = @(algebraicCandidate)residualSN_0(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:));
            soln = fsolve(residualfn, initialGuess_0, options);
            initialGuess_0 = soln;
            TEMP_CONSTANTS(:,128) = soln(1);
            TEMP_CONSTANTS(:,129) = soln(2);
            TEMP_CONSTANTS(:,130) = soln(3);
            ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC);
        end
    end
end

function resid = residualSN_0(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES)
    CONSTANTS(:,128) = algebraicCandidate(1);
    CONSTANTS(:,129) = algebraicCandidate(2);
    CONSTANTS(:,130) = algebraicCandidate(3);
    resid(1) = CONSTANTS(:,130) - ( CONSTANTS(:,128).*CONSTANTS(:,129))./CONSTANTS(:,29);
    resid(2) = CONSTANTS(:,128) - (CONSTANTS(:,16) - CONSTANTS(:,130));
    resid(3) = CONSTANTS(:,129) - (CONSTANTS(:,17) - CONSTANTS(:,130));
end

% Functions required for solving differential algebraic equation
function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_1(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN)
    ALGEBRAIC = ALGEBRAIC_IN;
    CONSTANTS = CONSTANTS_IN;
    STATES = STATES_IN;
    global initialGuess_1;
    if (length(initialGuess_1) ~= 17), initialGuess_1 = [0.1,0.1,0.1,0.988,5.5258e-05,3.8182,0.0588,0.227,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.0083];, end
    options = optimset('Display', 'off', 'TolX', 1E-6);
    if length(VOI) == 1
        residualfn = @(algebraicCandidate)residualSN_1(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES);
        soln = fsolve(residualfn, initialGuess_1, options);
        initialGuess_1 = soln;
        ALGEBRAIC(:,63) = soln(1);
        ALGEBRAIC(:,64) = soln(2);
        ALGEBRAIC(:,65) = soln(3);
        ALGEBRAIC(:,66) = soln(4);
        ALGEBRAIC(:,67) = soln(5);
        ALGEBRAIC(:,68) = soln(6);
        ALGEBRAIC(:,69) = soln(7);
        ALGEBRAIC(:,70) = soln(8);
        ALGEBRAIC(:,71) = soln(9);
        ALGEBRAIC(:,72) = soln(10);
        ALGEBRAIC(:,73) = soln(11);
        ALGEBRAIC(:,74) = soln(12);
        ALGEBRAIC(:,75) = soln(13);
        ALGEBRAIC(:,76) = soln(14);
        ALGEBRAIC(:,77) = soln(15);
        ALGEBRAIC(:,78) = soln(16);
        ALGEBRAIC(:,79) = soln(17);
    else
        SET_ALGEBRAIC(:,63) = logical(1);
        SET_ALGEBRAIC(:,64) = logical(1);
        SET_ALGEBRAIC(:,65) = logical(1);
        SET_ALGEBRAIC(:,66) = logical(1);
        SET_ALGEBRAIC(:,67) = logical(1);
        SET_ALGEBRAIC(:,68) = logical(1);
        SET_ALGEBRAIC(:,69) = logical(1);
        SET_ALGEBRAIC(:,70) = logical(1);
        SET_ALGEBRAIC(:,71) = logical(1);
        SET_ALGEBRAIC(:,72) = logical(1);
        SET_ALGEBRAIC(:,73) = logical(1);
        SET_ALGEBRAIC(:,74) = logical(1);
        SET_ALGEBRAIC(:,75) = logical(1);
        SET_ALGEBRAIC(:,76) = logical(1);
        SET_ALGEBRAIC(:,77) = logical(1);
        SET_ALGEBRAIC(:,78) = logical(1);
        SET_ALGEBRAIC(:,79) = logical(1);
        for i=1:length(VOI)
            residualfn = @(algebraicCandidate)residualSN_1(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:));
            soln = fsolve(residualfn, initialGuess_1, options);
            initialGuess_1 = soln;
            TEMP_ALGEBRAIC(:,63) = soln(1);
            TEMP_ALGEBRAIC(:,64) = soln(2);
            TEMP_ALGEBRAIC(:,65) = soln(3);
            TEMP_ALGEBRAIC(:,66) = soln(4);
            TEMP_ALGEBRAIC(:,67) = soln(5);
            TEMP_ALGEBRAIC(:,68) = soln(6);
            TEMP_ALGEBRAIC(:,69) = soln(7);
            TEMP_ALGEBRAIC(:,70) = soln(8);
            TEMP_ALGEBRAIC(:,71) = soln(9);
            TEMP_ALGEBRAIC(:,72) = soln(10);
            TEMP_ALGEBRAIC(:,73) = soln(11);
            TEMP_ALGEBRAIC(:,74) = soln(12);
            TEMP_ALGEBRAIC(:,75) = soln(13);
            TEMP_ALGEBRAIC(:,76) = soln(14);
            TEMP_ALGEBRAIC(:,77) = soln(15);
            TEMP_ALGEBRAIC(:,78) = soln(16);
            TEMP_ALGEBRAIC(:,79) = soln(17);
            ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC);
        end
    end
end

function resid = residualSN_1(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES)
    ALGEBRAIC(:,63) = algebraicCandidate(1);
    ALGEBRAIC(:,64) = algebraicCandidate(2);
    ALGEBRAIC(:,65) = algebraicCandidate(3);
    ALGEBRAIC(:,66) = algebraicCandidate(4);
    ALGEBRAIC(:,67) = algebraicCandidate(5);
    ALGEBRAIC(:,68) = algebraicCandidate(6);
    ALGEBRAIC(:,69) = algebraicCandidate(7);
    ALGEBRAIC(:,70) = algebraicCandidate(8);
    ALGEBRAIC(:,71) = algebraicCandidate(9);
    ALGEBRAIC(:,72) = algebraicCandidate(10);
    ALGEBRAIC(:,73) = algebraicCandidate(11);
    ALGEBRAIC(:,74) = algebraicCandidate(12);
    ALGEBRAIC(:,75) = algebraicCandidate(13);
    ALGEBRAIC(:,76) = algebraicCandidate(14);
    ALGEBRAIC(:,77) = algebraicCandidate(15);
    ALGEBRAIC(:,78) = algebraicCandidate(16);
    ALGEBRAIC(:,79) = algebraicCandidate(17);
    resid(1) = ALGEBRAIC(:,63) - ( ALGEBRAIC(:,66).*ALGEBRAIC(:,67))./CONSTANTS(:,4);
    resid(2) = ALGEBRAIC(:,64) - ( ALGEBRAIC(:,63).*ALGEBRAIC(:,68))./CONSTANTS(:,5);
    resid(3) = ALGEBRAIC(:,65) - ( ALGEBRAIC(:,67).*ALGEBRAIC(:,68))./CONSTANTS(:,6);
    resid(4) = ALGEBRAIC(:,66) - ((CONSTANTS(:,1) - ALGEBRAIC(:,63)) - ALGEBRAIC(:,64));
    resid(5) = ALGEBRAIC(:,67) - (((STATES(:,1) - ALGEBRAIC(:,63)) - ALGEBRAIC(:,64)) - ALGEBRAIC(:,65));
    resid(6) = ALGEBRAIC(:,68) - ((CONSTANTS(:,3) - ALGEBRAIC(:,64)) - ALGEBRAIC(:,65));
    resid(7) = ALGEBRAIC(:,71) - ( CONSTANTS(:,32).*CONSTANTS(:,36))./(CONSTANTS(:,36)+ALGEBRAIC(:,69)+ALGEBRAIC(:,79));
    resid(8) = ALGEBRAIC(:,72) -  (ALGEBRAIC(:,69)./CONSTANTS(:,35)).*ALGEBRAIC(:,69).*(1.00000+ALGEBRAIC(:,71)./CONSTANTS(:,36));
    resid(9) = ALGEBRAIC(:,73) -  ALGEBRAIC(:,69).*(1.00000+ALGEBRAIC(:,71)./CONSTANTS(:,36));
    resid(10) = ALGEBRAIC(:,74) -  (ALGEBRAIC(:,79)./CONSTANTS(:,35)).*ALGEBRAIC(:,79).*(1.00000+ALGEBRAIC(:,71)./CONSTANTS(:,36));
    resid(11) = ALGEBRAIC(:,75) -  ALGEBRAIC(:,79).*(1.00000+ALGEBRAIC(:,71)./CONSTANTS(:,36));
    resid(12) = ALGEBRAIC(:,76) -  (CONSTANTS(:,33)./ALGEBRAIC(:,70)).*ALGEBRAIC(:,72);
    resid(13) = ALGEBRAIC(:,77) -  (CONSTANTS(:,33)./ALGEBRAIC(:,70)).*ALGEBRAIC(:,74);
    resid(14) = ALGEBRAIC(:,78) - (( CONSTANTS(:,33).*CONSTANTS(:,34))./CONSTANTS(:,35)+( CONSTANTS(:,33).*ALGEBRAIC(:,70))./CONSTANTS(:,35)+power(ALGEBRAIC(:,70), 2.00000)./CONSTANTS(:,35));
    resid(15) = ALGEBRAIC(:,70) - ((STATES(:,7) - (ALGEBRAIC(:,76)+ 2.00000.*ALGEBRAIC(:,72)+ 2.00000.*ALGEBRAIC(:,73))) - (ALGEBRAIC(:,77)+ 2.00000.*ALGEBRAIC(:,74)+ 2.00000.*ALGEBRAIC(:,75)));
    resid(16) = 0.00000 - ( 2.00000.*CONSTANTS(:,30).*power(ALGEBRAIC(:,70), 2.00000) -  ALGEBRAIC(:,69).*(1.00000+ALGEBRAIC(:,71)./CONSTANTS(:,36)).*( ALGEBRAIC(:,78).*ALGEBRAIC(:,69)+power(ALGEBRAIC(:,70), 2.00000)));
    resid(17) = 0.00000 - ( 2.00000.*CONSTANTS(:,31).*power(ALGEBRAIC(:,70), 2.00000) -  ALGEBRAIC(:,79).*(1.00000+ALGEBRAIC(:,71)./CONSTANTS(:,36)).*( ALGEBRAIC(:,78).*ALGEBRAIC(:,79)+power(ALGEBRAIC(:,70), 2.00000)));
end

% Functions required for solving differential algebraic equation
function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_2(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN)
    ALGEBRAIC = ALGEBRAIC_IN;
    CONSTANTS = CONSTANTS_IN;
    STATES = STATES_IN;
    global initialGuess_2;
    if (length(initialGuess_2) ~= 3), initialGuess_2 = [0.0224,0.0471,0.1];, end
    options = optimset('Display', 'off', 'TolX', 1E-6);
    if length(VOI) == 1
        residualfn = @(algebraicCandidate)residualSN_2(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES);
        soln = fsolve(residualfn, initialGuess_2, options);
        initialGuess_2 = soln;
        ALGEBRAIC(:,43) = soln(1);
        ALGEBRAIC(:,44) = soln(2);
        ALGEBRAIC(:,45) = soln(3);
    else
        SET_ALGEBRAIC(:,43) = logical(1);
        SET_ALGEBRAIC(:,44) = logical(1);
        SET_ALGEBRAIC(:,45) = logical(1);
        for i=1:length(VOI)
            residualfn = @(algebraicCandidate)residualSN_2(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:));
            soln = fsolve(residualfn, initialGuess_2, options);
            initialGuess_2 = soln;
            TEMP_ALGEBRAIC(:,43) = soln(1);
            TEMP_ALGEBRAIC(:,44) = soln(2);
            TEMP_ALGEBRAIC(:,45) = soln(3);
            ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC);
        end
    end
end

function resid = residualSN_2(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES)
    ALGEBRAIC(:,43) = algebraicCandidate(1);
    ALGEBRAIC(:,44) = algebraicCandidate(2);
    ALGEBRAIC(:,45) = algebraicCandidate(3);
    resid(1) = ALGEBRAIC(:,45) - ( ALGEBRAIC(:,43).*ALGEBRAIC(:,44))./CONSTANTS(:,27);
    resid(2) = ALGEBRAIC(:,43) - (STATES(:,4) - ALGEBRAIC(:,45));
    resid(3) = ALGEBRAIC(:,44) - (CONSTANTS(:,14) - ALGEBRAIC(:,45));
end

% Functions required for solving differential algebraic equation
function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_3(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN)
    ALGEBRAIC = ALGEBRAIC_IN;
    CONSTANTS = CONSTANTS_IN;
    STATES = STATES_IN;
    global initialGuess_3;
    if (length(initialGuess_3) ~= 2), initialGuess_3 = [0,0.1];, end
    options = optimset('Display', 'off', 'TolX', 1E-6);
    if length(VOI) == 1
        residualfn = @(algebraicCandidate)residualSN_3(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES);
        soln = fsolve(residualfn, initialGuess_3, options);
        initialGuess_3 = soln;
        ALGEBRAIC(:,51) = soln(1);
        ALGEBRAIC(:,52) = soln(2);
    else
        SET_ALGEBRAIC(:,51) = logical(1);
        SET_ALGEBRAIC(:,52) = logical(1);
        for i=1:length(VOI)
            residualfn = @(algebraicCandidate)residualSN_3(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:));
            soln = fsolve(residualfn, initialGuess_3, options);
            initialGuess_3 = soln;
            TEMP_ALGEBRAIC(:,51) = soln(1);
            TEMP_ALGEBRAIC(:,52) = soln(2);
            ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC);
        end
    end
end

function resid = residualSN_3(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES)
    ALGEBRAIC(:,51) = algebraicCandidate(1);
    ALGEBRAIC(:,52) = algebraicCandidate(2);
    resid(1) = ALGEBRAIC(:,52) - ( ALGEBRAIC(:,51).*ALGEBRAIC(:,44))./CONSTANTS(:,28);
    resid(2) = ALGEBRAIC(:,51) - (CONSTANTS(:,18) - ALGEBRAIC(:,52));
end

% Functions required for solving differential algebraic equation
function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_4(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN)
    ALGEBRAIC = ALGEBRAIC_IN;
    CONSTANTS = CONSTANTS_IN;
    STATES = STATES_IN;
    global initialGuess_4;
    if (length(initialGuess_4) ~= 3), initialGuess_4 = [0.1,6.2734e-05,0.8375];, end
    options = optimset('Display', 'off', 'TolX', 1E-6);
    if length(VOI) == 1
        residualfn = @(algebraicCandidate)residualSN_4(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES);
        soln = fsolve(residualfn, initialGuess_4, options);
        initialGuess_4 = soln;
        ALGEBRAIC(:,56) = soln(1);
        ALGEBRAIC(:,57) = soln(2);
        ALGEBRAIC(:,58) = soln(3);
    else
        SET_ALGEBRAIC(:,56) = logical(1);
        SET_ALGEBRAIC(:,57) = logical(1);
        SET_ALGEBRAIC(:,58) = logical(1);
        for i=1:length(VOI)
            residualfn = @(algebraicCandidate)residualSN_4(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:));
            soln = fsolve(residualfn, initialGuess_4, options);
            initialGuess_4 = soln;
            TEMP_ALGEBRAIC(:,56) = soln(1);
            TEMP_ALGEBRAIC(:,57) = soln(2);
            TEMP_ALGEBRAIC(:,58) = soln(3);
            ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC);
        end
    end
end

function resid = residualSN_4(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES)
    ALGEBRAIC(:,56) = algebraicCandidate(1);
    ALGEBRAIC(:,57) = algebraicCandidate(2);
    ALGEBRAIC(:,58) = algebraicCandidate(3);
    resid(1) = ALGEBRAIC(:,56) - ( ALGEBRAIC(:,57).*ALGEBRAIC(:,58))./CONSTANTS(:,48);
    resid(2) = ALGEBRAIC(:,57) - (STATES(:,9) - ALGEBRAIC(:,56));
    resid(3) = ALGEBRAIC(:,58) - (CONSTANTS(:,38) - ALGEBRAIC(:,56));
end

% Compute result of a piecewise function
function x = piecewise(cases, default)
    set = [0];
    for i = 1:2:length(cases)
        if (length(cases{i+1}) == 1)
            x(cases{i} & ~set,:) = cases{i+1};
        else
            x(cases{i} & ~set,:) = cases{i+1}(cases{i} & ~set);
        end
        set = set | cases{i};
        if(set), break, end
    end
    if (length(default) == 1)
        x(~set,:) = default;
    else
        x(~set,:) = default(~set);
    end
end

% Pad out or shorten strings to a set length
function strout = strpad(strin)
    req_length = 160;
    insize = size(strin,2);
    if insize > req_length
        strout = strin(1:req_length);
    else
        strout = [strin, blanks(req_length - insize)];
    end
end