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 =120;
end
% There are a total of 41 entries in each of the rate and state variable arrays.
% There are a total of 137 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_CONSTANTS(:,1) = strpad('Bmax_CaM in component buffca (mM)');
    LEGEND_CONSTANTS(:,129) = strpad('Bmax_SLhighj in component buffca (mM)');
    LEGEND_CONSTANTS(:,134) = strpad('Bmax_SLhighsl in component buffca (mM)');
    LEGEND_CONSTANTS(:,130) = strpad('Bmax_SLlowj in component buffca (mM)');
    LEGEND_CONSTANTS(:,135) = strpad('Bmax_SLlowsl in component buffca (mM)');
    LEGEND_CONSTANTS(:,2) = strpad('Bmax_SR in component buffca (mM)');
    LEGEND_CONSTANTS(:,3) = strpad('Bmax_TnChigh in component buffca (mM)');
    LEGEND_CONSTANTS(:,4) = strpad('Bmax_TnClow in component buffca (mM)');
    LEGEND_CONSTANTS(:,5) = strpad('Bmax_myosin in component buffca (mM)');
    LEGEND_STATES(:,1) = strpad('CaM in component buffca (mM)');
    LEGEND_ALGEBRAIC(:,36) = strpad('J_CaB_cytosol in component buffca (mM_per_ms)');
    LEGEND_ALGEBRAIC(:,34) = strpad('J_CaB_jntion in component buffca (mM_per_ms)');
    LEGEND_ALGEBRAIC(:,35) = strpad('J_CaB_sl in component buffca (mM_per_ms)');
    LEGEND_STATES(:,2) = strpad('Myoc in component buffca (mM)');
    LEGEND_STATES(:,3) = strpad('Myom in component buffca (mM)');
    LEGEND_STATES(:,4) = strpad('SLH_jn in component buffca (mM)');
    LEGEND_STATES(:,5) = strpad('SLH_sl in component buffca (mM)');
    LEGEND_STATES(:,6) = strpad('SLL_jn in component buffca (mM)');
    LEGEND_STATES(:,7) = strpad('SLL_sl in component buffca (mM)');
    LEGEND_STATES(:,8) = strpad('SRB in component buffca (mM)');
    LEGEND_STATES(:,9) = strpad('TnCHc in component buffca (mM)');
    LEGEND_STATES(:,10) = strpad('TnCHm in component buffca (mM)');
    LEGEND_STATES(:,11) = strpad('TnCL in component buffca (mM)');
    LEGEND_CONSTANTS(:,6) = strpad('koff_cam in component buffca (mS_per_uF)');
    LEGEND_CONSTANTS(:,7) = strpad('koff_myoca in component buffca (mS_per_uF)');
    LEGEND_CONSTANTS(:,8) = strpad('koff_myomg in component buffca (mS_per_uF)');
    LEGEND_CONSTANTS(:,9) = strpad('koff_slh in component buffca (mS_per_uF)');
    LEGEND_CONSTANTS(:,10) = strpad('koff_sll in component buffca (mS_per_uF)');
    LEGEND_CONSTANTS(:,11) = strpad('koff_sr in component buffca (mS_per_uF)');
    LEGEND_CONSTANTS(:,12) = strpad('koff_tnchca in component buffca (mS_per_uF)');
    LEGEND_CONSTANTS(:,13) = strpad('koff_tnchmg in component buffca (mS_per_uF)');
    LEGEND_CONSTANTS(:,103) = strpad('koff_tncl in component buffca (mS_per_uF)');
    LEGEND_CONSTANTS(:,14) = strpad('kon_cam in component buffca (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,15) = strpad('kon_myoca in component buffca (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,16) = strpad('kon_myomg in component buffca (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,17) = strpad('kon_slh in component buffca (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,18) = strpad('kon_sll in component buffca (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,19) = strpad('kon_sr in component buffca (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,20) = strpad('kon_tnchca in component buffca (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,21) = strpad('kon_tnchmg in component buffca (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,22) = strpad('kon_tncl in component buffca (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,124) = strpad('Vjn in component geom (liter)');
    LEGEND_CONSTANTS(:,127) = strpad('Vmyo in component geom (liter)');
    LEGEND_CONSTANTS(:,131) = strpad('Vsl in component geom (liter)');
    LEGEND_STATES(:,12) = strpad('Ca_i in component calcium (mM)');
    LEGEND_VOI = strpad('time in component engine (ms)');
    LEGEND_CONSTANTS(:,23) = strpad('Mg_i in component ion (mM)');
    LEGEND_STATES(:,13) = strpad('Ca_jn in component calcium (mM)');
    LEGEND_STATES(:,14) = strpad('Ca_sl in component calcium (mM)');
    LEGEND_CONSTANTS(:,24) = strpad('ISO in component cell (dimensionless)');
    LEGEND_CONSTANTS(:,25) = strpad('Bmax_Na_jn in component buffna (mM)');
    LEGEND_CONSTANTS(:,26) = strpad('Bmax_Na_sl in component buffna (mM)');
    LEGEND_STATES(:,15) = strpad('NaB_jn in component buffna (mM)');
    LEGEND_STATES(:,16) = strpad('NaB_sl in component buffna (mM)');
    LEGEND_CONSTANTS(:,27) = strpad('koff_na in component buffna (mS_per_uF)');
    LEGEND_CONSTANTS(:,28) = strpad('kon_na in component buffna (per_mM_per_ms)');
    LEGEND_STATES(:,17) = strpad('Na_jn in component sodium (mM)');
    LEGEND_STATES(:,18) = strpad('Na_sl in component sodium (mM)');
    LEGEND_STATES(:,19) = strpad('Ca_sr in component calcium (mM)');
    LEGEND_STATES(:,20) = strpad('Csqn in component calcium (mM)');
    LEGEND_ALGEBRAIC(:,83) = strpad('ICa_tot_jn in component calcium (A_per_F)');
    LEGEND_ALGEBRAIC(:,90) = strpad('ICa_tot_sl in component calcium (A_per_F)');
    LEGEND_CONSTANTS(:,29) = strpad('JCa_slmyo in component geom (m3_per_s)');
    LEGEND_CONSTANTS(:,132) = strpad('Vsr in component geom (liter)');
    LEGEND_ALGEBRAIC(:,88) = strpad('J_serca in component ryr (mM_per_ms)');
    LEGEND_CONSTANTS(:,30) = strpad('C in component cell (farad)');
    LEGEND_CONSTANTS(:,31) = strpad('JCa_jnsl in component geom (m3_per_s)');
    LEGEND_CONSTANTS(:,32) = strpad('F in component phys (C_per_mol)');
    LEGEND_ALGEBRAIC(:,84) = strpad('J_SRCarel in component ryr (mM_per_ms)');
    LEGEND_ALGEBRAIC(:,86) = strpad('J_SR_leak in component ryr (mM_per_ms)');
    LEGEND_CONSTANTS(:,136) = strpad('Bmax_Csqn in component calcium_Csqn (mM)');
    LEGEND_CONSTANTS(:,33) = strpad('koff_csqn in component calcium_Csqn (mS_per_uF)');
    LEGEND_CONSTANTS(:,34) = strpad('kon_csqn in component calcium_Csqn (per_mM_per_ms)');
    LEGEND_ALGEBRAIC(:,82) = strpad('ICaB_jn in component icab (A_per_F)');
    LEGEND_ALGEBRAIC(:,38) = strpad('ICaL_Ca_jn in component ical (A_per_F)');
    LEGEND_ALGEBRAIC(:,65) = strpad('INaCa_jn in component inaca (A_per_F)');
    LEGEND_ALGEBRAIC(:,76) = strpad('IpCa_jn in component ipca (A_per_F)');
    LEGEND_ALGEBRAIC(:,87) = strpad('ICaB_sl in component icab (A_per_F)');
    LEGEND_ALGEBRAIC(:,40) = strpad('ICaL_Ca_sl in component ical (A_per_F)');
    LEGEND_ALGEBRAIC(:,69) = strpad('INaCa_sl in component inaca (A_per_F)');
    LEGEND_ALGEBRAIC(:,78) = strpad('IpCa_sl in component ipca (A_per_F)');
    LEGEND_CONSTANTS(:,35) = strpad('AF in component cell (dimensionless)');
    LEGEND_CONSTANTS(:,36) = strpad('RA in component cell (dimensionless)');
    LEGEND_CONSTANTS(:,37) = strpad('JNa_jnsl in component geom (m3_per_s)');
    LEGEND_CONSTANTS(:,38) = strpad('JNa_slmyo in component geom (m3_per_s)');
    LEGEND_CONSTANTS(:,104) = strpad('Vcell in component geom (liter)');
    LEGEND_CONSTANTS(:,39) = strpad('cell_length in component geom (um)');
    LEGEND_CONSTANTS(:,40) = strpad('cell_radius in component geom (um)');
    LEGEND_CONSTANTS(:,41) = strpad('pi in component geom (dimensionless)');
    LEGEND_ALGEBRAIC(:,89) = strpad('ICaB in component icab (A_per_F)');
    LEGEND_CONSTANTS(:,42) = strpad('gCaB in component icab (mS_per_uF)');
    LEGEND_CONSTANTS(:,43) = strpad('Fjn in component junc (dimensionless)');
    LEGEND_STATES(:,21) = strpad('V in component membrane (mV)');
    LEGEND_ALGEBRAIC(:,81) = strpad('ECa_jn in component nernst (mV)');
    LEGEND_CONSTANTS(:,128) = strpad('Fsl in component junc (dimensionless)');
    LEGEND_ALGEBRAIC(:,85) = strpad('ECa_sl in component nernst (mV)');
    LEGEND_ALGEBRAIC(:,51) = strpad('ICaL in component ical (A_per_F)');
    LEGEND_ALGEBRAIC(:,41) = strpad('ICaL_Ca in component ical (A_per_F)');
    LEGEND_ALGEBRAIC(:,43) = strpad('ICaL_K in component ical (A_per_F)');
    LEGEND_ALGEBRAIC(:,49) = strpad('ICaL_Na in component ical (A_per_F)');
    LEGEND_ALGEBRAIC(:,45) = strpad('ICaL_Na_jn in component ical (A_per_F)');
    LEGEND_ALGEBRAIC(:,47) = strpad('ICaL_Na_sl in component ical (A_per_F)');
    LEGEND_CONSTANTS(:,44) = strpad('Q10CaL in component ical (dimensionless)');
    LEGEND_STATES(:,22) = strpad('d in component ical (dimensionless)');
    LEGEND_STATES(:,23) = strpad('f in component ical (dimensionless)');
    LEGEND_STATES(:,24) = strpad('fCaB_jn in component ical (dimensionless)');
    LEGEND_STATES(:,25) = strpad('fCaB_sl in component ical (dimensionless)');
    LEGEND_CONSTANTS(:,45) = strpad('f_conducting in component ical (dimensionless)');
    LEGEND_CONSTANTS(:,46) = strpad('fcaCaMSL in component ical (dimensionless)');
    LEGEND_CONSTANTS(:,47) = strpad('fcaCaj in component ical (dimensionless)');
    LEGEND_ALGEBRAIC(:,37) = strpad('ibarca_jn in component ical (uA_per_cm2)');
    LEGEND_ALGEBRAIC(:,39) = strpad('ibarca_sl in component ical (uA_per_cm2)');
    LEGEND_ALGEBRAIC(:,42) = strpad('ibark in component ical (uA_per_cm2)');
    LEGEND_ALGEBRAIC(:,44) = strpad('ibarna_jn in component ical (uA_per_cm2)');
    LEGEND_ALGEBRAIC(:,46) = strpad('ibarna_sl in component ical (uA_per_cm2)');
    LEGEND_CONSTANTS(:,106) = strpad('pCa in component ical (cm_per_s)');
    LEGEND_CONSTANTS(:,48) = strpad('pCa_max in component ical (cm_per_s)');
    LEGEND_CONSTANTS(:,107) = strpad('pK in component ical (cm_per_s)');
    LEGEND_CONSTANTS(:,49) = strpad('pK_max in component ical (cm_per_s)');
    LEGEND_CONSTANTS(:,108) = strpad('pNa in component ical (cm_per_s)');
    LEGEND_CONSTANTS(:,50) = strpad('pNa_max in component ical (cm_per_s)');
    LEGEND_CONSTANTS(:,51) = strpad('Fjn_CaL in component junc (dimensionless)');
    LEGEND_CONSTANTS(:,133) = strpad('Q in component phys (dimensionless)');
    LEGEND_CONSTANTS(:,109) = strpad('Fsl_CaL in component junc (dimensionless)');
    LEGEND_ALGEBRAIC(:,1) = strpad('ical_d_inf in component ical_d (dimensionless)');
    LEGEND_ALGEBRAIC(:,15) = strpad('ical_d_tau in component ical_d (ms)');
    LEGEND_ALGEBRAIC(:,2) = strpad('ical_f_inf in component ical_f (dimensionless)');
    LEGEND_ALGEBRAIC(:,16) = strpad('ical_f_tau in component ical_f (ms)');
    LEGEND_CONSTANTS(:,52) = strpad('Ca_o in component ion (mM)');
    LEGEND_CONSTANTS(:,110) = strpad('FRT in component phys (per_mV)');
    LEGEND_CONSTANTS(:,53) = strpad('K_o in component ion (mM)');
    LEGEND_STATES(:,26) = strpad('K_i in component potassium (mM)');
    LEGEND_CONSTANTS(:,54) = strpad('Na_o in component ion (mM)');
    LEGEND_ALGEBRAIC(:,48) = strpad('IClB in component iclb (A_per_F)');
    LEGEND_CONSTANTS(:,55) = strpad('gClB in component iclb (mS_per_uF)');
    LEGEND_CONSTANTS(:,126) = strpad('ECl in component nernst (mV)');
    LEGEND_CONSTANTS(:,56) = strpad('GClCa in component iclca (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,53) = strpad('IClCa in component iclca (A_per_F)');
    LEGEND_ALGEBRAIC(:,50) = strpad('IClCa_jn in component iclca (A_per_F)');
    LEGEND_ALGEBRAIC(:,52) = strpad('IClCa_sl in component iclca (A_per_F)');
    LEGEND_CONSTANTS(:,57) = strpad('KdClCa in component iclca (mM)');
    LEGEND_ALGEBRAIC(:,96) = strpad('IK1 in component ik1 (A_per_F)');
    LEGEND_CONSTANTS(:,111) = strpad('gK1 in component ik1 (mS_per_uF)');
    LEGEND_CONSTANTS(:,58) = strpad('gK1_max in component ik1 (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,95) = strpad('ik1_inf in component ik1 (dimensionless)');
    LEGEND_ALGEBRAIC(:,92) = strpad('EK in component nernst (mV)');
    LEGEND_ALGEBRAIC(:,93) = strpad('ik1_inf_a in component ik1_ik1_inf (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,94) = strpad('ik1_inf_b in component ik1_ik1_inf (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,99) = strpad('IKp in component ikp (A_per_F)');
    LEGEND_ALGEBRAIC(:,97) = strpad('IKp_jn in component ikp (A_per_F)');
    LEGEND_ALGEBRAIC(:,98) = strpad('IKp_sl in component ikp (A_per_F)');
    LEGEND_CONSTANTS(:,59) = strpad('gKp in component ikp (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,54) = strpad('kp in component ikp (dimensionless)');
    LEGEND_ALGEBRAIC(:,100) = strpad('IKr in component ikr (A_per_F)');
    LEGEND_CONSTANTS(:,112) = strpad('gKr in component ikr (mS_per_uF)');
    LEGEND_CONSTANTS(:,60) = strpad('gKr_max in component ikr (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,55) = strpad('rr in component ikr (dimensionless)');
    LEGEND_STATES(:,27) = strpad('xr in component ikr (dimensionless)');
    LEGEND_ALGEBRAIC(:,3) = strpad('ikr_xr_inf in component ikr_xr (dimensionless)');
    LEGEND_ALGEBRAIC(:,17) = strpad('ikr_xr_tau in component ikr_xr (ms)');
    LEGEND_ALGEBRAIC(:,56) = strpad('EKs in component iks (mV)');
    LEGEND_ALGEBRAIC(:,59) = strpad('IKs in component iks (A_per_F)');
    LEGEND_ALGEBRAIC(:,57) = strpad('IKs_jn in component iks (A_per_F)');
    LEGEND_ALGEBRAIC(:,58) = strpad('IKs_sl in component iks (A_per_F)');
    LEGEND_CONSTANTS(:,113) = strpad('gKs_jn in component iks (mS_per_uF)');
    LEGEND_CONSTANTS(:,61) = strpad('gKs_max in component iks (mS_per_uF)');
    LEGEND_CONSTANTS(:,114) = strpad('gKs_sl in component iks (mS_per_uF)');
    LEGEND_CONSTANTS(:,62) = strpad('pNaK in component iks (dimensionless)');
    LEGEND_STATES(:,28) = strpad('xs in component iks (dimensionless)');
    LEGEND_STATES(:,29) = strpad('Na_i in component sodium (mM)');
    LEGEND_ALGEBRAIC(:,4) = strpad('iks_xs_inf in component iks_xs (dimensionless)');
    LEGEND_ALGEBRAIC(:,18) = strpad('iks_xs_tau in component iks_xs (ms)');
    LEGEND_ALGEBRAIC(:,101) = strpad('IKur in component ikur (A_per_F)');
    LEGEND_CONSTANTS(:,115) = strpad('gKur in component ikur (mS_per_uF)');
    LEGEND_CONSTANTS(:,63) = strpad('gKur_max in component ikur (mS_per_uF)');
    LEGEND_STATES(:,30) = strpad('ikur_r in component ikur (dimensionless)');
    LEGEND_STATES(:,31) = strpad('s in component ikur (dimensionless)');
    LEGEND_ALGEBRAIC(:,5) = strpad('ikur_r_inf in component ikur_ikur_r (dimensionless)');
    LEGEND_ALGEBRAIC(:,19) = strpad('ikur_r_tau in component ikur_ikur_r (ms)');
    LEGEND_ALGEBRAIC(:,6) = strpad('ikur_s_inf in component ikur_s (dimensionless)');
    LEGEND_ALGEBRAIC(:,20) = strpad('ikur_s_tau in component ikur_s (ms)');
    LEGEND_ALGEBRAIC(:,109) = strpad('INa in component ina (A_per_F)');
    LEGEND_ALGEBRAIC(:,104) = strpad('INa_jn in component ina (A_per_F)');
    LEGEND_ALGEBRAIC(:,108) = strpad('INa_sl in component ina (A_per_F)');
    LEGEND_CONSTANTS(:,116) = strpad('gNa in component ina (mS_per_uF)');
    LEGEND_CONSTANTS(:,64) = strpad('gNa_max in component ina (mS_per_uF)');
    LEGEND_STATES(:,32) = strpad('h in component ina (dimensionless)');
    LEGEND_STATES(:,33) = strpad('j in component ina (dimensionless)');
    LEGEND_STATES(:,34) = strpad('m in component ina (dimensionless)');
    LEGEND_ALGEBRAIC(:,103) = strpad('ENa_jn in component nernst (mV)');
    LEGEND_ALGEBRAIC(:,107) = strpad('ENa_sl in component nernst (mV)');
    LEGEND_ALGEBRAIC(:,28) = strpad('ina_h_inf in component ina_h (dimensionless)');
    LEGEND_ALGEBRAIC(:,31) = strpad('ina_h_tau in component ina_h (ms)');
    LEGEND_ALGEBRAIC(:,29) = strpad('ina_j_inf in component ina_j (dimensionless)');
    LEGEND_ALGEBRAIC(:,32) = strpad('ina_j_tau in component ina_j (ms)');
    LEGEND_ALGEBRAIC(:,9) = strpad('ina_m_inf in component ina_m (dimensionless)');
    LEGEND_ALGEBRAIC(:,23) = strpad('ina_m_tau in component ina_m (ms)');
    LEGEND_ALGEBRAIC(:,7) = strpad('ina_h_a in component ina_h (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,21) = strpad('ina_h_b in component ina_h (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,8) = strpad('ina_j_a in component ina_j (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,22) = strpad('ina_j_b in component ina_j (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,111) = strpad('INaB in component inab (A_per_F)');
    LEGEND_ALGEBRAIC(:,105) = strpad('INaB_jn in component inab (A_per_F)');
    LEGEND_ALGEBRAIC(:,110) = strpad('INaB_sl in component inab (A_per_F)');
    LEGEND_CONSTANTS(:,65) = strpad('gNaB in component inab (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,71) = strpad('INaCa in component inaca (A_per_F)');
    LEGEND_CONSTANTS(:,117) = strpad('IbarNCX in component inaca (A_per_F)');
    LEGEND_CONSTANTS(:,66) = strpad('IbarNCX_max in component inaca (A_per_F)');
    LEGEND_ALGEBRAIC(:,60) = strpad('Ka_jn in component inaca (dimensionless)');
    LEGEND_ALGEBRAIC(:,61) = strpad('Ka_sl in component inaca (dimensionless)');
    LEGEND_CONSTANTS(:,67) = strpad('Kdact in component inaca (mM)');
    LEGEND_CONSTANTS(:,68) = strpad('KmCai in component inaca (mM)');
    LEGEND_CONSTANTS(:,69) = strpad('KmCao in component inaca (mM)');
    LEGEND_CONSTANTS(:,70) = strpad('KmNai in component inaca (mM)');
    LEGEND_CONSTANTS(:,71) = strpad('KmNao in component inaca (mM)');
    LEGEND_CONSTANTS(:,72) = strpad('Q10NCX in component inaca (dimensionless)');
    LEGEND_CONSTANTS(:,73) = strpad('ksat in component inaca (dimensionless)');
    LEGEND_CONSTANTS(:,74) = strpad('nu in component inaca (dimensionless)');
    LEGEND_ALGEBRAIC(:,62) = strpad('s1_jn in component inaca_INaCa_jn (mol4_per_m12)');
    LEGEND_ALGEBRAIC(:,63) = strpad('s2_jn in component inaca_INaCa_jn (mol4_per_m12)');
    LEGEND_ALGEBRAIC(:,64) = strpad('s3_jn in component inaca_INaCa_jn (mol4_per_m12)');
    LEGEND_ALGEBRAIC(:,66) = strpad('s1_sl in component inaca_INaCa_sl (mol4_per_m12)');
    LEGEND_ALGEBRAIC(:,67) = strpad('s2_sl in component inaca_INaCa_sl (mol4_per_m12)');
    LEGEND_ALGEBRAIC(:,68) = strpad('s3_sl in component inaca_INaCa_sl (mol4_per_m12)');
    LEGEND_ALGEBRAIC(:,74) = strpad('INaK in component inak (A_per_F)');
    LEGEND_ALGEBRAIC(:,72) = strpad('INaK_jn in component inak (A_per_F)');
    LEGEND_ALGEBRAIC(:,73) = strpad('INaK_sl in component inak (A_per_F)');
    LEGEND_CONSTANTS(:,75) = strpad('IbarNaK in component inak (A_per_F)');
    LEGEND_CONSTANTS(:,76) = strpad('KmKo in component inak (mM)');
    LEGEND_CONSTANTS(:,118) = strpad('KmNaip in component inak (mM)');
    LEGEND_ALGEBRAIC(:,70) = strpad('fnak in component inak (dimensionless)');
    LEGEND_CONSTANTS(:,119) = strpad('sigma in component inak (dimensionless)');
    LEGEND_ALGEBRAIC(:,113) = strpad('INaL in component inal (A_per_F)');
    LEGEND_ALGEBRAIC(:,106) = strpad('INaL_jn in component inal (A_per_F)');
    LEGEND_ALGEBRAIC(:,112) = strpad('INaL_sl in component inal (A_per_F)');
    LEGEND_CONSTANTS(:,120) = strpad('gNaL in component inal (mS_per_uF)');
    LEGEND_CONSTANTS(:,77) = strpad('gNaL_max in component inal (mS_per_uF)');
    LEGEND_STATES(:,35) = strpad('hl in component inal (dimensionless)');
    LEGEND_STATES(:,36) = strpad('ml in component inal (dimensionless)');
    LEGEND_ALGEBRAIC(:,10) = strpad('inal_hl_inf in component inal_hl (dimensionless)');
    LEGEND_CONSTANTS(:,78) = strpad('inal_hl_tau in component inal_hl (ms)');
    LEGEND_ALGEBRAIC(:,11) = strpad('inal_ml_a in component inal_ml (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,24) = strpad('inal_ml_b in component inal_ml (mS_per_uF)');
    LEGEND_CONSTANTS(:,79) = strpad('Cl_i in component ion (mM)');
    LEGEND_CONSTANTS(:,80) = strpad('Cl_o in component ion (mM)');
    LEGEND_CONSTANTS(:,81) = strpad('IbarSLCaP in component ipca (A_per_F)');
    LEGEND_ALGEBRAIC(:,80) = strpad('IpCa in component ipca (A_per_F)');
    LEGEND_CONSTANTS(:,82) = strpad('KmPCa in component ipca (mM)');
    LEGEND_CONSTANTS(:,83) = strpad('Q10SLCaP in component ipca (dimensionless)');
    LEGEND_ALGEBRAIC(:,75) = strpad('ipca_IpCa_jn_a in component ipca_IpCa_jn (dimensionless)');
    LEGEND_CONSTANTS(:,105) = strpad('ipca_IpCa_jn_b in component ipca_IpCa_jn (dimensionless)');
    LEGEND_ALGEBRAIC(:,77) = strpad('ipca_IpCa_sl_a in component ipca_IpCa_sl (dimensionless)');
    LEGEND_CONSTANTS(:,125) = strpad('ipca_IpCa_sl_b in component ipca_IpCa_sl (dimensionless)');
    LEGEND_ALGEBRAIC(:,102) = strpad('Ito in component ito (A_per_F)');
    LEGEND_CONSTANTS(:,121) = strpad('gto in component ito (mS_per_uF)');
    LEGEND_CONSTANTS(:,84) = strpad('gto_max in component ito (mS_per_uF)');
    LEGEND_STATES(:,37) = strpad('x in component ito (dimensionless)');
    LEGEND_STATES(:,38) = strpad('y in component ito (dimensionless)');
    LEGEND_ALGEBRAIC(:,12) = strpad('ito_x_inf in component ito_x (dimensionless)');
    LEGEND_ALGEBRAIC(:,25) = strpad('ito_x_tau in component ito_x (ms)');
    LEGEND_ALGEBRAIC(:,13) = strpad('ito_y_inf in component ito_y (dimensionless)');
    LEGEND_ALGEBRAIC(:,26) = strpad('ito_y_tau in component ito_y (ms)');
    LEGEND_ALGEBRAIC(:,91) = strpad('ICa_tot in component membrane (A_per_F)');
    LEGEND_ALGEBRAIC(:,79) = strpad('ICl_tot in component membrane (A_per_F)');
    LEGEND_ALGEBRAIC(:,117) = strpad('INa_tot in component membrane (A_per_F)');
    LEGEND_ALGEBRAIC(:,118) = strpad('I_tot in component membrane (A_per_F)');
    LEGEND_ALGEBRAIC(:,115) = strpad('INa_tot_jn in component sodium (A_per_F)');
    LEGEND_ALGEBRAIC(:,116) = strpad('INa_tot_sl in component sodium (A_per_F)');
    LEGEND_ALGEBRAIC(:,114) = strpad('IK_tot in component potassium (A_per_F)');
    LEGEND_ALGEBRAIC(:,120) = strpad('I_stim in component stimulus (A_per_F)');
    LEGEND_CONSTANTS(:,85) = strpad('R in component phys (mJ_per_mol_per_K)');
    LEGEND_CONSTANTS(:,86) = strpad('T in component phys (kelvin)');
    LEGEND_CONSTANTS(:,87) = strpad('J_SR_leak_max in component ryr (mS_per_uF)');
    LEGEND_CONSTANTS(:,122) = strpad('Kmf in component ryr (mM)');
    LEGEND_CONSTANTS(:,88) = strpad('Kmr in component ryr (mM)');
    LEGEND_CONSTANTS(:,89) = strpad('MaxSR in component ryr (dimensionless)');
    LEGEND_CONSTANTS(:,90) = strpad('MinSR in component ryr (dimensionless)');
    LEGEND_CONSTANTS(:,91) = strpad('Q10SRCaP in component ryr (dimensionless)');
    LEGEND_ALGEBRAIC(:,14) = strpad('RI in component ryr (dimensionless)');
    LEGEND_CONSTANTS(:,92) = strpad('Vmax_SRCaP in component ryr (mM_per_ms)');
    LEGEND_CONSTANTS(:,93) = strpad('ec50SR in component ryr (mM)');
    LEGEND_CONSTANTS(:,94) = strpad('hillSRCaP in component ryr (dimensionless)');
    LEGEND_STATES(:,39) = strpad('i in component ryr (dimensionless)');
    LEGEND_ALGEBRAIC(:,27) = strpad('kCaSR in component ryr (dimensionless)');
    LEGEND_CONSTANTS(:,95) = strpad('kiCa in component ryr (per_mM_per_ms)');
    LEGEND_ALGEBRAIC(:,30) = strpad('kiSRCa in component ryr (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,96) = strpad('kim in component ryr (mS_per_uF)');
    LEGEND_CONSTANTS(:,123) = strpad('koCa in component ryr (per_mM2_per_ms)');
    LEGEND_ALGEBRAIC(:,33) = strpad('koSRCa in component ryr (per_mM2_per_ms)');
    LEGEND_CONSTANTS(:,97) = strpad('kom in component ryr (mS_per_uF)');
    LEGEND_CONSTANTS(:,98) = strpad('ks in component ryr (mS_per_uF)');
    LEGEND_STATES(:,40) = strpad('o in component ryr (dimensionless)');
    LEGEND_STATES(:,41) = strpad('ryr_r in component ryr (dimensionless)');
    LEGEND_CONSTANTS(:,99) = strpad('amplitude in component stimulus (A_per_F)');
    LEGEND_ALGEBRAIC(:,119) = strpad('pace in component stimulus (dimensionless)');
    LEGEND_CONSTANTS(:,100) = strpad('stim_duration in component stimulus (ms)');
    LEGEND_CONSTANTS(:,101) = strpad('stim_offset in component stimulus (ms)');
    LEGEND_CONSTANTS(:,102) = strpad('stim_period in component stimulus (ms)');
    LEGEND_RATES(:,1) = strpad('d/dt CaM in component buffca (mM)');
    LEGEND_RATES(:,11) = strpad('d/dt TnCL in component buffca (mM)');
    LEGEND_RATES(:,9) = strpad('d/dt TnCHc in component buffca (mM)');
    LEGEND_RATES(:,10) = strpad('d/dt TnCHm in component buffca (mM)');
    LEGEND_RATES(:,2) = strpad('d/dt Myoc in component buffca (mM)');
    LEGEND_RATES(:,3) = strpad('d/dt Myom in component buffca (mM)');
    LEGEND_RATES(:,8) = strpad('d/dt SRB in component buffca (mM)');
    LEGEND_RATES(:,6) = strpad('d/dt SLL_jn in component buffca (mM)');
    LEGEND_RATES(:,4) = strpad('d/dt SLH_jn in component buffca (mM)');
    LEGEND_RATES(:,7) = strpad('d/dt SLL_sl in component buffca (mM)');
    LEGEND_RATES(:,5) = strpad('d/dt SLH_sl in component buffca (mM)');
    LEGEND_RATES(:,15) = strpad('d/dt NaB_jn in component buffna (mM)');
    LEGEND_RATES(:,16) = strpad('d/dt NaB_sl in component buffna (mM)');
    LEGEND_RATES(:,12) = strpad('d/dt Ca_i in component calcium (mM)');
    LEGEND_RATES(:,13) = strpad('d/dt Ca_jn in component calcium (mM)');
    LEGEND_RATES(:,14) = strpad('d/dt Ca_sl in component calcium (mM)');
    LEGEND_RATES(:,19) = strpad('d/dt Ca_sr in component calcium (mM)');
    LEGEND_RATES(:,20) = strpad('d/dt Csqn in component calcium (mM)');
    LEGEND_RATES(:,22) = strpad('d/dt d in component ical (dimensionless)');
    LEGEND_RATES(:,23) = strpad('d/dt f in component ical (dimensionless)');
    LEGEND_RATES(:,24) = strpad('d/dt fCaB_jn in component ical (dimensionless)');
    LEGEND_RATES(:,25) = strpad('d/dt fCaB_sl in component ical (dimensionless)');
    LEGEND_RATES(:,27) = strpad('d/dt xr in component ikr (dimensionless)');
    LEGEND_RATES(:,28) = strpad('d/dt xs in component iks (dimensionless)');
    LEGEND_RATES(:,30) = strpad('d/dt ikur_r in component ikur (dimensionless)');
    LEGEND_RATES(:,31) = strpad('d/dt s in component ikur (dimensionless)');
    LEGEND_RATES(:,32) = strpad('d/dt h in component ina (dimensionless)');
    LEGEND_RATES(:,33) = strpad('d/dt j in component ina (dimensionless)');
    LEGEND_RATES(:,34) = strpad('d/dt m in component ina (dimensionless)');
    LEGEND_RATES(:,35) = strpad('d/dt hl in component inal (dimensionless)');
    LEGEND_RATES(:,36) = strpad('d/dt ml in component inal (dimensionless)');
    LEGEND_RATES(:,37) = strpad('d/dt x in component ito (dimensionless)');
    LEGEND_RATES(:,38) = strpad('d/dt y in component ito (dimensionless)');
    LEGEND_RATES(:,21) = strpad('d/dt V in component membrane (mV)');
    LEGEND_RATES(:,26) = strpad('d/dt K_i in component potassium (mM)');
    LEGEND_RATES(:,39) = strpad('d/dt i in component ryr (dimensionless)');
    LEGEND_RATES(:,40) = strpad('d/dt o in component ryr (dimensionless)');
    LEGEND_RATES(:,41) = strpad('d/dt ryr_r in component ryr (dimensionless)');
    LEGEND_RATES(:,29) = strpad('d/dt Na_i in component sodium (mM)');
    LEGEND_RATES(:,17) = strpad('d/dt Na_jn in component sodium (mM)');
    LEGEND_RATES(:,18) = strpad('d/dt Na_sl in component sodium (mM)');
    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) = 0.024;
    CONSTANTS(:,2) = 0.0171;
    CONSTANTS(:,3) = 0.14;
    CONSTANTS(:,4) = 0.07;
    CONSTANTS(:,5) = 0.14;
    STATES(:,1) =  7.02128101897185673e-04;
    STATES(:,2) =  3.94923428392655786e-03;
    STATES(:,3) =  1.35538532457244482e-01;
    STATES(:,4) =  1.03674364292988680e-01;
    STATES(:,5) =  1.90759804527589089e-01;
    STATES(:,6) =  1.35640688636079511e-02;
    STATES(:,7) =  2.14063418881809235e-02;
    STATES(:,8) =  4.45327242854324807e-03;
    STATES(:,9) =  1.27856586024588575e-01;
    STATES(:,10) =  5.69999505293381902e-03;
    STATES(:,11) =  1.83143535034222225e-02;
    CONSTANTS(:,6) = 0.238;
    CONSTANTS(:,7) = 0.00046;
    CONSTANTS(:,8) = 5.7e-05;
    CONSTANTS(:,9) = 0.03;
    CONSTANTS(:,10) = 1.3;
    CONSTANTS(:,11) = 0.06;
    CONSTANTS(:,12) = 3.2e-05;
    CONSTANTS(:,13) = 0.00333;
    CONSTANTS(:,14) = 34.0;
    CONSTANTS(:,15) = 13.8;
    CONSTANTS(:,16) = 0.0157;
    CONSTANTS(:,17) = 100.0;
    CONSTANTS(:,18) = 100.0;
    CONSTANTS(:,19) = 100.0;
    CONSTANTS(:,20) = 2.37;
    CONSTANTS(:,21) = 0.003;
    CONSTANTS(:,22) = 32.7;
    STATES(:,12) =  2.10808768153058460e-04;
    CONSTANTS(:,23) = 1.0;
    STATES(:,13) =  3.25814677291117296e-04;
    STATES(:,14) =  2.33018340557575125e-04;
    CONSTANTS(:,24) = 0.0;
    CONSTANTS(:,25) = 7.561;
    CONSTANTS(:,26) = 1.65;
    STATES(:,15) =  3.61396062660070427e+00;
    STATES(:,16) =  7.88607791910409195e-01;
    CONSTANTS(:,27) = 0.001;
    CONSTANTS(:,28) = 0.0001;
    STATES(:,17) =  9.15153381546177336e+00;
    STATES(:,18) =  9.15182798281732346e+00;
    STATES(:,19) =  5.02305826642838293e-01;
    STATES(:,20) =  1.13337536953687845e+00;
    CONSTANTS(:,29) =  3.72425607984805052e-12;
    CONSTANTS(:,30) = 1.1e-10;
    CONSTANTS(:,31) =  8.24130542277896849e-13;
    CONSTANTS(:,32) = 96485.0;
    CONSTANTS(:,33) = 65.0;
    CONSTANTS(:,34) = 100.0;
    CONSTANTS(:,35) = 0.0;
    CONSTANTS(:,36) = 0.0;
    CONSTANTS(:,37) =  1.83127823220607955e-14;
    CONSTANTS(:,38) =  1.63862792221979433e-12;
    CONSTANTS(:,39) = 100.0;
    CONSTANTS(:,40) = 10.25;
    CONSTANTS(:,41) =  3.14159265358979312e+00;
    CONSTANTS(:,42) =  6.06430000000000033e-04;
    CONSTANTS(:,43) = 0.11;
    STATES(:,21) = -7.34336366728778671e+01;
    CONSTANTS(:,44) = 1.8;
    STATES(:,22) =  2.16850216379767157e-05;
    STATES(:,23) =  9.98384427312367095e-01;
    STATES(:,24) =  4.49572164109603364e-02;
    STATES(:,25) =  3.28512098597005947e-02;
    CONSTANTS(:,45) = 1.0;
    CONSTANTS(:,46) = 0.0;
    CONSTANTS(:,47) = 0.0;
    CONSTANTS(:,48) = 0.00027;
    CONSTANTS(:,49) = 1.35e-07;
    CONSTANTS(:,50) = 7.5e-09;
    CONSTANTS(:,51) = 0.9;
    CONSTANTS(:,52) = 1.8;
    CONSTANTS(:,53) = 5.4;
    STATES(:,26) = 120.0;
    CONSTANTS(:,54) = 140.0;
    CONSTANTS(:,55) = 0.009;
    CONSTANTS(:,56) = 0.0548;
    CONSTANTS(:,57) = 0.1;
    CONSTANTS(:,58) = 0.0525;
    CONSTANTS(:,59) = 0.002;
    CONSTANTS(:,60) = 0.035;
    STATES(:,27) =  1.31290096227093382e-03;
    CONSTANTS(:,61) = 0.0035;
    CONSTANTS(:,62) = 0.01833;
    STATES(:,28) =  7.49436760722081534e-03;
    STATES(:,29) =  9.15199678386256998e+00;
    CONSTANTS(:,63) = 0.045;
    STATES(:,30) =  3.93548562883350357e-04;
    STATES(:,31) =  9.58234428284286399e-01;
    CONSTANTS(:,64) = 23.0;
    STATES(:,32) =  3.15482710277587786e-01;
    STATES(:,33) =  2.48034071360795916e-01;
    STATES(:,34) =  1.89326933812916480e-02;
    CONSTANTS(:,65) = 0.000597;
    CONSTANTS(:,66) = 3.15;
    CONSTANTS(:,67) = 0.000384;
    CONSTANTS(:,68) = 0.00359;
    CONSTANTS(:,69) = 1.3;
    CONSTANTS(:,70) = 12.29;
    CONSTANTS(:,71) = 87.5;
    CONSTANTS(:,72) = 1.57;
    CONSTANTS(:,73) = 0.27;
    CONSTANTS(:,74) = 0.35;
    CONSTANTS(:,75) = 1.26;
    CONSTANTS(:,76) = 1.5;
    CONSTANTS(:,77) = 0.0025;
    STATES(:,35) =  3.79829335413739144e-02;
    STATES(:,36) =  1.01974216400706526e-02;
    CONSTANTS(:,78) = 600.0;
    CONSTANTS(:,79) = 15.0;
    CONSTANTS(:,80) = 150.0;
    CONSTANTS(:,81) = 0.0471;
    CONSTANTS(:,82) = 0.0005;
    CONSTANTS(:,83) = 2.35;
    CONSTANTS(:,84) = 0.165;
    STATES(:,37) =  1.37939236359928058e-03;
    STATES(:,38) =  9.45874848392074696e-01;
    CONSTANTS(:,85) = 8314.0;
    CONSTANTS(:,86) = 310.0;
    CONSTANTS(:,87) = 5.348e-06;
    CONSTANTS(:,88) = 1.7;
    CONSTANTS(:,89) = 15.0;
    CONSTANTS(:,90) = 1.0;
    CONSTANTS(:,91) = 2.6;
    CONSTANTS(:,92) = 0.0053114;
    CONSTANTS(:,93) = 0.45;
    CONSTANTS(:,94) = 1.787;
    STATES(:,39) =  5.01323282772066123e-07;
    CONSTANTS(:,95) = 0.5;
    CONSTANTS(:,96) = 0.005;
    CONSTANTS(:,97) = 0.06;
    CONSTANTS(:,98) = 25.0;
    STATES(:,40) =  2.01567245823636694e-06;
    STATES(:,41) =  8.00819151705148946e-01;
    CONSTANTS(:,99) = -12.5;
    CONSTANTS(:,100) = 5.0;
    CONSTANTS(:,101) = 50.0;
    CONSTANTS(:,102) = 1000.0;
    CONSTANTS(:,103) =  (1.00000+ 0.500000.*CONSTANTS(:,24)).*0.0196000;
    CONSTANTS(:,104) =  ( ( CONSTANTS(:,41).*power(CONSTANTS(:,40), 2.00000)).*CONSTANTS(:,39)).*1.00000e-15;
    CONSTANTS(:,105) = power( CONSTANTS(:,82).*1.00000, 1.60000);
    CONSTANTS(:,106) =  ( (1.00000+ 0.500000.*CONSTANTS(:,24)).*(1.00000 -  0.500000.*CONSTANTS(:,35))).*CONSTANTS(:,48);
    CONSTANTS(:,107) =  ( (1.00000+ 0.500000.*CONSTANTS(:,24)).*(1.00000 -  0.500000.*CONSTANTS(:,35))).*CONSTANTS(:,49);
    CONSTANTS(:,108) =  ( (1.00000+ 0.500000.*CONSTANTS(:,24)).*(1.00000 -  0.500000.*CONSTANTS(:,35))).*CONSTANTS(:,50);
    CONSTANTS(:,109) = 1.00000 - CONSTANTS(:,51);
    CONSTANTS(:,110) = (CONSTANTS(:,32)./CONSTANTS(:,85))./CONSTANTS(:,86);
    CONSTANTS(:,111) =  ( (1.00000+CONSTANTS(:,35)).*power((CONSTANTS(:,53)./5.40000), 1.0 ./ 2)).*CONSTANTS(:,58);
    CONSTANTS(:,112) =  CONSTANTS(:,60).*power((CONSTANTS(:,53)./5.40000), 1.0 ./ 2);
    CONSTANTS(:,113) =  ((1.00000+CONSTANTS(:,35))+ 2.00000.*CONSTANTS(:,24)).*CONSTANTS(:,61);
    CONSTANTS(:,114) =  ((1.00000+CONSTANTS(:,35))+ 2.00000.*CONSTANTS(:,24)).*CONSTANTS(:,61);
    CONSTANTS(:,115) =  ( ( (1.00000 -  0.500000.*CONSTANTS(:,35)).*(1.00000+ 2.00000.*CONSTANTS(:,24))).*(1.00000+ 0.200000.*CONSTANTS(:,36))).*CONSTANTS(:,63);
    CONSTANTS(:,116) =  CONSTANTS(:,64).*(1.00000 -  0.100000.*CONSTANTS(:,35));
    CONSTANTS(:,117) =  (1.00000+ 0.400000.*CONSTANTS(:,35)).*CONSTANTS(:,66);
    CONSTANTS(:,118) =  11.0000.*(1.00000 -  0.250000.*CONSTANTS(:,24));
    CONSTANTS(:,119) = (exp(CONSTANTS(:,54)./67.3000) - 1.00000)./7.00000;
    CONSTANTS(:,120) =  CONSTANTS(:,77).*CONSTANTS(:,35);
    CONSTANTS(:,121) =  (1.00000 -  0.700000.*CONSTANTS(:,35)).*CONSTANTS(:,84);
    CONSTANTS(:,122) =  (2.50000 -  1.25000.*CONSTANTS(:,24)).*0.000246000;
    CONSTANTS(:,123) =  ((10.0000+ 20.0000.*CONSTANTS(:,35))+ ( 10.0000.*CONSTANTS(:,24)).*(1.00000 - CONSTANTS(:,35))).*1.00000;
    CONSTANTS(:,136) = 0.00000;
    CONSTANTS(:,124) =  ( 0.0539000.*0.0100000).*CONSTANTS(:,104);
    CONSTANTS(:,125) = power( CONSTANTS(:,82).*1.00000, 1.60000);
    CONSTANTS(:,126) =  (1.00000./CONSTANTS(:,110)).*log(CONSTANTS(:,79)./CONSTANTS(:,80));
    CONSTANTS(:,127) =  0.650000.*CONSTANTS(:,104);
    CONSTANTS(:,128) = 1.00000 - CONSTANTS(:,43);
    CONSTANTS(:,129) =  (( 0.00165000.*CONSTANTS(:,127))./CONSTANTS(:,124)).*0.100000;
    CONSTANTS(:,130) =  (( 0.00460000.*CONSTANTS(:,127))./CONSTANTS(:,124)).*0.100000;
    CONSTANTS(:,131) =  0.0200000.*CONSTANTS(:,104);
    CONSTANTS(:,132) =  0.0350000.*CONSTANTS(:,104);
    CONSTANTS(:,133) = (CONSTANTS(:,86) - 310.000)./10.0000;
    CONSTANTS(:,134) =  (CONSTANTS(:,127)./CONSTANTS(:,131)).*0.0134000;
    CONSTANTS(:,135) =  (CONSTANTS(:,127)./CONSTANTS(:,131)).*0.0374000;
    CONSTANTS(:,136) =  (CONSTANTS(:,127)./CONSTANTS(:,132)).*0.140000;
    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
    RATES(:,26) = CONSTANTS(:,136);
    RATES(:,1) =  ( CONSTANTS(:,14).*STATES(:,12)).*(CONSTANTS(:,1) - STATES(:,1)) -  CONSTANTS(:,6).*STATES(:,1);
    RATES(:,24) =  ( ( 1.70000.*STATES(:,13)).*(1.00000 - STATES(:,24)) -  0.0119000.*STATES(:,24)).*1.00000;
    RATES(:,25) =  ( ( 1.70000.*STATES(:,14)).*(1.00000 - STATES(:,25)) -  0.0119000.*STATES(:,25)).*1.00000;
    RATES(:,29) =  (CONSTANTS(:,38)./CONSTANTS(:,127)).*(STATES(:,18) - STATES(:,29));
    RATES(:,2) =  ( CONSTANTS(:,15).*STATES(:,12)).*((CONSTANTS(:,5) - STATES(:,2)) - STATES(:,3)) -  CONSTANTS(:,7).*STATES(:,2);
    ALGEBRAIC(:,10) = 1.00000./(1.00000+exp((STATES(:,21)+91.0000)./6.10000));
    RATES(:,35) = (ALGEBRAIC(:,10) - STATES(:,35))./CONSTANTS(:,78);
    RATES(:,3) =  ( CONSTANTS(:,16).*CONSTANTS(:,23)).*((CONSTANTS(:,5) - STATES(:,2)) - STATES(:,3)) -  CONSTANTS(:,8).*STATES(:,3);
    ALGEBRAIC(:,1) = 1.00000./(1.00000+exp( - ((STATES(:,21)+ 3.00000.*CONSTANTS(:,24))+9.00000)./6.00000));
    ALGEBRAIC(:,15) = ( ALGEBRAIC(:,1).*(1.00000 - exp( - ((STATES(:,21)+ 3.00000.*CONSTANTS(:,24))+9.00000)./6.00000)))./( 0.0350000.*((STATES(:,21)+ 3.00000.*CONSTANTS(:,24))+9.00000));
    RATES(:,22) = (ALGEBRAIC(:,1) - STATES(:,22))./ALGEBRAIC(:,15);
    ALGEBRAIC(:,2) = 1.00000./(1.00000+exp(((STATES(:,21)+ 3.00000.*CONSTANTS(:,24))+30.0000)./7.00000))+0.200000./(1.00000+exp(((50.0000 - STATES(:,21)) -  3.00000.*CONSTANTS(:,24))./20.0000));
    ALGEBRAIC(:,16) = 1.00000./( 0.0197000.*exp( - power( 0.0337000.*((STATES(:,21)+ 3.00000.*CONSTANTS(:,24))+25.0000), 2.00000))+0.0200000);
    RATES(:,23) = (ALGEBRAIC(:,2) - STATES(:,23))./ALGEBRAIC(:,16);
    ALGEBRAIC(:,3) = 1.00000./(1.00000+exp( - (STATES(:,21)+10.0000)./5.00000));
    ALGEBRAIC(:,17) = ( ((550.000)./(1.00000+exp(( - 22.0000 - STATES(:,21))./9.00000))).*6.00000)./(1.00000+exp((STATES(:,21)+11.0000)./9.00000))+230.000./(1.00000+exp((STATES(:,21)+40.0000)./20.0000));
    RATES(:,27) = (ALGEBRAIC(:,3) - STATES(:,27))./ALGEBRAIC(:,17);
    ALGEBRAIC(:,4) = 1.00000./(1.00000+exp( - ((STATES(:,21)+ 40.0000.*CONSTANTS(:,24))+3.80000)./14.2500));
    ALGEBRAIC(:,18) = 990.100./(1.00000+exp( - ((STATES(:,21)+ 40.0000.*CONSTANTS(:,24))+2.43600)./14.1200));
    RATES(:,28) = (ALGEBRAIC(:,4) - STATES(:,28))./ALGEBRAIC(:,18);
    ALGEBRAIC(:,5) = 1.00000./(1.00000+exp((STATES(:,21)+6.00000)./ - 8.60000));
    ALGEBRAIC(:,19) = 9.00000./(1.00000+exp((STATES(:,21)+5.00000)./12.0000))+0.500000;
    RATES(:,30) = (ALGEBRAIC(:,5) - STATES(:,30))./ALGEBRAIC(:,19);
    ALGEBRAIC(:,6) = 1.00000./(1.00000+exp((STATES(:,21)+7.50000)./10.0000));
    ALGEBRAIC(:,20) = 590.000./(1.00000+exp((STATES(:,21)+60.0000)./10.0000))+3050.00;
    RATES(:,31) = (ALGEBRAIC(:,6) - STATES(:,31))./ALGEBRAIC(:,20);
    ALGEBRAIC(:,9) = 1.00000./power(1.00000+exp( - (56.8600+STATES(:,21))./9.03000), 2.00000);
    ALGEBRAIC(:,23) =  0.129200.*exp( - power((STATES(:,21)+45.7900)./15.5400, 2.00000))+ 0.0648700.*exp( - power((STATES(:,21) - 4.82300)./51.1200, 2.00000));
    RATES(:,34) = (ALGEBRAIC(:,9) - STATES(:,34))./ALGEBRAIC(:,23);
    ALGEBRAIC(:,11) = ( 0.320000.*(STATES(:,21)+47.1300))./(1.00000 - exp(  - 0.100000.*(STATES(:,21)+47.1300)));
    ALGEBRAIC(:,24) =  0.0800000.*exp( - STATES(:,21)./11.0000);
    RATES(:,36) =  ALGEBRAIC(:,11).*(1.00000 - STATES(:,36)) -  ALGEBRAIC(:,24).*STATES(:,36);
    ALGEBRAIC(:,12) = 1.00000./(1.00000+exp( - (STATES(:,21)+1.00000)./11.0000));
    ALGEBRAIC(:,25) =  3.50000.*exp( - power(STATES(:,21)./30.0000, 2.00000))+1.50000;
    RATES(:,37) = (ALGEBRAIC(:,12) - STATES(:,37))./ALGEBRAIC(:,25);
    ALGEBRAIC(:,13) = 1.00000./(1.00000+exp((STATES(:,21)+40.5000)./11.5000));
    ALGEBRAIC(:,26) =  25.6350.*exp( - power((STATES(:,21)+52.4500)./15.8827, 2.00000))+24.1400;
    RATES(:,38) = (ALGEBRAIC(:,13) - STATES(:,38))./ALGEBRAIC(:,26);
    RATES(:,4) =  ( CONSTANTS(:,17).*STATES(:,13)).*(CONSTANTS(:,129) - STATES(:,4)) -  CONSTANTS(:,9).*STATES(:,4);
    RATES(:,5) =  ( CONSTANTS(:,17).*STATES(:,14)).*(CONSTANTS(:,134) - STATES(:,5)) -  CONSTANTS(:,9).*STATES(:,5);
    ALGEBRAIC(:,28) = 1.00000./power(1.00000+exp((STATES(:,21)+71.5500)./7.43000), 2.00000);
    ALGEBRAIC(:,7) =  piecewise({STATES(:,21)>= - 40.0000, 0.00000 },  0.0570000.*exp( - (STATES(:,21)+80.0000)./6.80000)).*1.00000;
    ALGEBRAIC(:,21) =  piecewise({STATES(:,21)>= - 40.0000, 0.770000./( 0.130000.*(1.00000+exp( - (STATES(:,21)+10.6600)./11.1000))) },  2.70000.*exp( 0.0790000.*STATES(:,21))+ ( 3.10000.*power(10.0000, 5.00000)).*exp( 0.348500.*STATES(:,21))).*1.00000;
    ALGEBRAIC(:,31) = 1.00000./(ALGEBRAIC(:,7)+ALGEBRAIC(:,21));
    RATES(:,32) = (ALGEBRAIC(:,28) - STATES(:,32))./ALGEBRAIC(:,31);
    ALGEBRAIC(:,29) = 1.00000./power(1.00000+exp((STATES(:,21)+71.5500)./7.43000), 2.00000);
    ALGEBRAIC(:,8) =  piecewise({STATES(:,21)>= - 40.0000, 0.00000 },  (( ( (  - 2.54280.*power(10.0000, 4.00000)).*exp( 0.244400.*STATES(:,21)) -  ( 6.94800.*power(10.0000,  - 6.00000)).*exp(  - 0.0439100.*STATES(:,21))).*(STATES(:,21)+37.7800))./(1.00000+exp( 0.311000.*(STATES(:,21)+79.2300)))).*1.00000).*1.00000;
    ALGEBRAIC(:,22) =  piecewise({STATES(:,21)>= - 40.0000, ( 0.600000.*exp( 0.0570000.*STATES(:,21)))./(1.00000+exp(  - 0.100000.*(STATES(:,21)+32.0000))) }, ( 0.0242400.*exp(  - 0.0105200.*STATES(:,21)))./(1.00000+exp(  - 0.137800.*(STATES(:,21)+40.1400)))).*1.00000;
    ALGEBRAIC(:,32) = 1.00000./(ALGEBRAIC(:,8)+ALGEBRAIC(:,22));
    RATES(:,33) = (ALGEBRAIC(:,29) - STATES(:,33))./ALGEBRAIC(:,32);
    ALGEBRAIC(:,14) = ((1.00000 - STATES(:,41)) - STATES(:,40)) - STATES(:,39);
    ALGEBRAIC(:,27) = CONSTANTS(:,89) - (CONSTANTS(:,89) - CONSTANTS(:,90))./(1.00000+power(CONSTANTS(:,93)./STATES(:,19), 2.50000));
    ALGEBRAIC(:,30) =  CONSTANTS(:,95).*ALGEBRAIC(:,27);
    ALGEBRAIC(:,33) = CONSTANTS(:,123)./ALGEBRAIC(:,27);
    RATES(:,39) = ( ( ALGEBRAIC(:,30).*STATES(:,13)).*STATES(:,40) -  CONSTANTS(:,96).*STATES(:,39)) - ( CONSTANTS(:,97).*STATES(:,39) -  ( ALGEBRAIC(:,33).*power(STATES(:,13), 2.00000)).*ALGEBRAIC(:,14));
    RATES(:,40) = ( ( ALGEBRAIC(:,33).*power(STATES(:,13), 2.00000)).*STATES(:,41) -  CONSTANTS(:,97).*STATES(:,40)) - ( ( ALGEBRAIC(:,30).*STATES(:,13)).*STATES(:,40) -  CONSTANTS(:,96).*STATES(:,39));
    RATES(:,41) = ( CONSTANTS(:,96).*ALGEBRAIC(:,14) -  ( ALGEBRAIC(:,30).*STATES(:,13)).*STATES(:,41)) - ( ( ALGEBRAIC(:,33).*power(STATES(:,13), 2.00000)).*STATES(:,41) -  CONSTANTS(:,97).*STATES(:,40));
    RATES(:,6) =  ( CONSTANTS(:,18).*STATES(:,13)).*(CONSTANTS(:,130) - STATES(:,6)) -  CONSTANTS(:,10).*STATES(:,6);
    RATES(:,7) =  ( CONSTANTS(:,18).*STATES(:,14)).*(CONSTANTS(:,135) - STATES(:,7)) -  CONSTANTS(:,10).*STATES(:,7);
    RATES(:,8) =  ( CONSTANTS(:,19).*STATES(:,12)).*(CONSTANTS(:,2) - STATES(:,8)) -  CONSTANTS(:,11).*STATES(:,8);
    RATES(:,9) =  ( CONSTANTS(:,20).*STATES(:,12)).*((CONSTANTS(:,3) - STATES(:,9)) - STATES(:,10)) -  CONSTANTS(:,12).*STATES(:,9);
    RATES(:,10) =  ( CONSTANTS(:,21).*CONSTANTS(:,23)).*((CONSTANTS(:,3) - STATES(:,9)) - STATES(:,10)) -  CONSTANTS(:,13).*STATES(:,10);
    RATES(:,11) =  ( CONSTANTS(:,22).*STATES(:,12)).*(CONSTANTS(:,4) - STATES(:,11)) -  CONSTANTS(:,103).*STATES(:,11);
    RATES(:,15) =  ( CONSTANTS(:,28).*STATES(:,17)).*(CONSTANTS(:,25) - STATES(:,15)) -  CONSTANTS(:,27).*STATES(:,15);
    RATES(:,16) =  ( CONSTANTS(:,28).*STATES(:,18)).*(CONSTANTS(:,26) - STATES(:,16)) -  CONSTANTS(:,27).*STATES(:,16);
    RATES(:,20) =  ( CONSTANTS(:,34).*STATES(:,19)).*(CONSTANTS(:,136) - STATES(:,20)) -  CONSTANTS(:,33).*STATES(:,20);
    ALGEBRAIC(:,34) = RATES(:,6)+RATES(:,4);
    ALGEBRAIC(:,81) =  ((1.00000./CONSTANTS(:,110))./2.00000).*log(CONSTANTS(:,52)./STATES(:,13));
    ALGEBRAIC(:,82) =  ( CONSTANTS(:,43).*CONSTANTS(:,42)).*(STATES(:,21) - ALGEBRAIC(:,81));
    ALGEBRAIC(:,37) =  (( ( ( CONSTANTS(:,106).*4.00000).*( ( STATES(:,21).*CONSTANTS(:,32)).*CONSTANTS(:,110))).*( ( 0.341000.*STATES(:,13)).*exp( ( 2.00000.*STATES(:,21)).*CONSTANTS(:,110)) -  0.341000.*CONSTANTS(:,52)))./(exp( ( 2.00000.*STATES(:,21)).*CONSTANTS(:,110)) - 1.00000)).*CONSTANTS(:,45);
    ALGEBRAIC(:,38) =  ( ( ( ( ( CONSTANTS(:,51).*ALGEBRAIC(:,37)).*STATES(:,22)).*STATES(:,23)).*((1.00000 - STATES(:,24))+CONSTANTS(:,47))).*power(CONSTANTS(:,44), CONSTANTS(:,133))).*0.450000;
    ALGEBRAIC(:,60) = 1.00000./(1.00000+power(CONSTANTS(:,67)./STATES(:,13), 2.00000));
    ALGEBRAIC(:,62) =  ( exp( ( CONSTANTS(:,74).*STATES(:,21)).*CONSTANTS(:,110)).*power(STATES(:,17), 3.00000)).*CONSTANTS(:,52);
    ALGEBRAIC(:,63) =  ( exp( ( (CONSTANTS(:,74) - 1.00000).*STATES(:,21)).*CONSTANTS(:,110)).*power(CONSTANTS(:,54), 3.00000)).*STATES(:,13);
    ALGEBRAIC(:,64) = ((( ( CONSTANTS(:,68).*power(CONSTANTS(:,54), 3.00000)).*(1.00000+power(STATES(:,17)./CONSTANTS(:,70), 3.00000))+ ( power(CONSTANTS(:,71), 3.00000).*STATES(:,13)).*(1.00000+STATES(:,13)./CONSTANTS(:,68)))+ CONSTANTS(:,69).*power(STATES(:,17), 3.00000))+ power(STATES(:,17), 3.00000).*CONSTANTS(:,52))+ power(CONSTANTS(:,54), 3.00000).*STATES(:,13);
    ALGEBRAIC(:,65) = (( ( ( ( CONSTANTS(:,43).*CONSTANTS(:,117)).*power(CONSTANTS(:,72), CONSTANTS(:,133))).*ALGEBRAIC(:,60)).*(ALGEBRAIC(:,62) - ALGEBRAIC(:,63)))./ALGEBRAIC(:,64))./(1.00000+ CONSTANTS(:,73).*exp( ( (CONSTANTS(:,74) - 1.00000).*STATES(:,21)).*CONSTANTS(:,110)));
    ALGEBRAIC(:,75) = power( STATES(:,13).*1.00000, 1.60000);
    ALGEBRAIC(:,76) = ( ( ( CONSTANTS(:,43).*power(CONSTANTS(:,83), CONSTANTS(:,133))).*CONSTANTS(:,81)).*ALGEBRAIC(:,75))./(ALGEBRAIC(:,75)+CONSTANTS(:,105));
    ALGEBRAIC(:,83) = ((ALGEBRAIC(:,38)+ALGEBRAIC(:,82))+ALGEBRAIC(:,76)) -  2.00000.*ALGEBRAIC(:,65);
    ALGEBRAIC(:,84) =  ( CONSTANTS(:,98).*STATES(:,40)).*(STATES(:,19) - STATES(:,13));
    ALGEBRAIC(:,86) =  ( (1.00000+ 0.250000.*CONSTANTS(:,35)).*(STATES(:,19) - STATES(:,13))).*CONSTANTS(:,87);
    RATES(:,13) = ((((  - ALGEBRAIC(:,83).*CONSTANTS(:,30))./( ( CONSTANTS(:,124).*2.00000).*CONSTANTS(:,32))+ (CONSTANTS(:,31)./CONSTANTS(:,124)).*(STATES(:,14) - STATES(:,13))) - ALGEBRAIC(:,34))+( ALGEBRAIC(:,84).*CONSTANTS(:,132))./CONSTANTS(:,124))+( ALGEBRAIC(:,86).*CONSTANTS(:,127))./CONSTANTS(:,124);
    ALGEBRAIC(:,36) = (((((RATES(:,11)+RATES(:,9))+RATES(:,10))+RATES(:,1))+RATES(:,2))+RATES(:,3))+RATES(:,8);
    ALGEBRAIC(:,88) = ( ( power(CONSTANTS(:,91), CONSTANTS(:,133)).*CONSTANTS(:,92)).*(power(STATES(:,12)./CONSTANTS(:,122), CONSTANTS(:,94)) - power(STATES(:,19)./CONSTANTS(:,88), CONSTANTS(:,94))))./((1.00000+power(STATES(:,12)./CONSTANTS(:,122), CONSTANTS(:,94)))+power(STATES(:,19)./CONSTANTS(:,88), CONSTANTS(:,94)));
    RATES(:,12) = ((  - ALGEBRAIC(:,88).*CONSTANTS(:,132))./CONSTANTS(:,127) - ALGEBRAIC(:,36))+ (CONSTANTS(:,29)./CONSTANTS(:,127)).*(STATES(:,14) - STATES(:,12));
    ALGEBRAIC(:,35) = RATES(:,7)+RATES(:,5);
    ALGEBRAIC(:,85) =  ((1.00000./CONSTANTS(:,110))./2.00000).*log(CONSTANTS(:,52)./STATES(:,14));
    ALGEBRAIC(:,87) =  ( CONSTANTS(:,128).*CONSTANTS(:,42)).*(STATES(:,21) - ALGEBRAIC(:,85));
    ALGEBRAIC(:,39) =  (( ( ( CONSTANTS(:,106).*4.00000).*( ( STATES(:,21).*CONSTANTS(:,32)).*CONSTANTS(:,110))).*( ( 0.341000.*STATES(:,14)).*exp( ( 2.00000.*STATES(:,21)).*CONSTANTS(:,110)) -  0.341000.*CONSTANTS(:,52)))./(exp( ( 2.00000.*STATES(:,21)).*CONSTANTS(:,110)) - 1.00000)).*CONSTANTS(:,45);
    ALGEBRAIC(:,40) =  ( ( ( ( ( CONSTANTS(:,109).*ALGEBRAIC(:,39)).*STATES(:,22)).*STATES(:,23)).*((1.00000 - STATES(:,25))+CONSTANTS(:,46))).*power(CONSTANTS(:,44), CONSTANTS(:,133))).*0.450000;
    ALGEBRAIC(:,61) = 1.00000./(1.00000+power(CONSTANTS(:,67)./STATES(:,14), 2.00000));
    ALGEBRAIC(:,66) =  ( exp( ( CONSTANTS(:,74).*STATES(:,21)).*CONSTANTS(:,110)).*power(STATES(:,18), 3.00000)).*CONSTANTS(:,52);
    ALGEBRAIC(:,67) =  ( exp( ( (CONSTANTS(:,74) - 1.00000).*STATES(:,21)).*CONSTANTS(:,110)).*power(CONSTANTS(:,54), 3.00000)).*STATES(:,14);
    ALGEBRAIC(:,68) = ((( ( CONSTANTS(:,68).*power(CONSTANTS(:,54), 3.00000)).*(1.00000+power(STATES(:,18)./CONSTANTS(:,70), 3.00000))+ ( power(CONSTANTS(:,71), 3.00000).*STATES(:,14)).*(1.00000+STATES(:,14)./CONSTANTS(:,68)))+ CONSTANTS(:,69).*power(STATES(:,18), 3.00000))+ power(STATES(:,18), 3.00000).*CONSTANTS(:,52))+ power(CONSTANTS(:,54), 3.00000).*STATES(:,14);
    ALGEBRAIC(:,69) = (( ( ( ( CONSTANTS(:,128).*CONSTANTS(:,117)).*power(CONSTANTS(:,72), CONSTANTS(:,133))).*ALGEBRAIC(:,61)).*(ALGEBRAIC(:,66) - ALGEBRAIC(:,67)))./ALGEBRAIC(:,68))./(1.00000+ CONSTANTS(:,73).*exp( ( (CONSTANTS(:,74) - 1.00000).*STATES(:,21)).*CONSTANTS(:,110)));
    ALGEBRAIC(:,77) = power( STATES(:,14).*1.00000, 1.60000);
    ALGEBRAIC(:,78) = ( ( ( CONSTANTS(:,128).*power(CONSTANTS(:,83), CONSTANTS(:,133))).*CONSTANTS(:,81)).*ALGEBRAIC(:,77))./(ALGEBRAIC(:,77)+CONSTANTS(:,125));
    ALGEBRAIC(:,90) = ((ALGEBRAIC(:,40)+ALGEBRAIC(:,87))+ALGEBRAIC(:,78)) -  2.00000.*ALGEBRAIC(:,69);
    RATES(:,14) = (((  - ALGEBRAIC(:,90).*CONSTANTS(:,30))./( ( CONSTANTS(:,131).*2.00000).*CONSTANTS(:,32))+ (CONSTANTS(:,31)./CONSTANTS(:,131)).*(STATES(:,13) - STATES(:,14)))+ (CONSTANTS(:,29)./CONSTANTS(:,131)).*(STATES(:,12) - STATES(:,14))) - ALGEBRAIC(:,35);
    RATES(:,19) = (ALGEBRAIC(:,88) - (( ALGEBRAIC(:,86).*CONSTANTS(:,127))./CONSTANTS(:,132)+ALGEBRAIC(:,84))) - RATES(:,20);
    ALGEBRAIC(:,44) =  (( ( ( ( CONSTANTS(:,108).*STATES(:,21)).*CONSTANTS(:,32)).*CONSTANTS(:,110)).*( ( 0.750000.*STATES(:,17)).*exp( STATES(:,21).*CONSTANTS(:,110)) -  0.750000.*CONSTANTS(:,54)))./(exp( STATES(:,21).*CONSTANTS(:,110)) - 1.00000)).*CONSTANTS(:,45);
    ALGEBRAIC(:,45) =  ( ( ( ( ( CONSTANTS(:,51).*ALGEBRAIC(:,44)).*STATES(:,22)).*STATES(:,23)).*((1.00000 - STATES(:,24))+CONSTANTS(:,47))).*power(CONSTANTS(:,44), CONSTANTS(:,133))).*0.450000;
    ALGEBRAIC(:,103) =  (1.00000./CONSTANTS(:,110)).*log(CONSTANTS(:,54)./STATES(:,17));
    ALGEBRAIC(:,104) =  ( ( ( ( CONSTANTS(:,43).*CONSTANTS(:,116)).*power(STATES(:,34), 3.00000)).*STATES(:,32)).*STATES(:,33)).*(STATES(:,21) - ALGEBRAIC(:,103));
    ALGEBRAIC(:,105) =  ( CONSTANTS(:,43).*CONSTANTS(:,65)).*(STATES(:,21) - ALGEBRAIC(:,103));
    ALGEBRAIC(:,70) = 1.00000./((1.00000+ 0.124500.*exp( (  - 0.100000.*STATES(:,21)).*CONSTANTS(:,110)))+ ( 0.0365000.*CONSTANTS(:,119)).*exp(  - STATES(:,21).*CONSTANTS(:,110)));
    ALGEBRAIC(:,72) = (( ( ( CONSTANTS(:,43).*CONSTANTS(:,75)).*ALGEBRAIC(:,70)).*CONSTANTS(:,53))./(1.00000+power(CONSTANTS(:,118)./STATES(:,17), 4.00000)))./(CONSTANTS(:,53)+CONSTANTS(:,76));
    ALGEBRAIC(:,106) =  ( ( ( CONSTANTS(:,43).*CONSTANTS(:,120)).*power(STATES(:,36), 3.00000)).*STATES(:,35)).*(STATES(:,21) - ALGEBRAIC(:,103));
    ALGEBRAIC(:,115) = ((((ALGEBRAIC(:,104)+ALGEBRAIC(:,105))+ 3.00000.*ALGEBRAIC(:,65))+ 3.00000.*ALGEBRAIC(:,72))+ALGEBRAIC(:,45))+ALGEBRAIC(:,106);
    RATES(:,17) = ((  - ALGEBRAIC(:,115).*CONSTANTS(:,30))./( CONSTANTS(:,124).*CONSTANTS(:,32))+ (CONSTANTS(:,37)./CONSTANTS(:,124)).*(STATES(:,18) - STATES(:,17))) - RATES(:,15);
    ALGEBRAIC(:,46) =  (( ( ( ( CONSTANTS(:,108).*STATES(:,21)).*CONSTANTS(:,32)).*CONSTANTS(:,110)).*( ( 0.750000.*STATES(:,18)).*exp( STATES(:,21).*CONSTANTS(:,110)) -  0.750000.*CONSTANTS(:,54)))./(exp( STATES(:,21).*CONSTANTS(:,110)) - 1.00000)).*CONSTANTS(:,45);
    ALGEBRAIC(:,47) =  ( ( ( ( ( CONSTANTS(:,109).*ALGEBRAIC(:,46)).*STATES(:,22)).*STATES(:,23)).*((1.00000 - STATES(:,25))+CONSTANTS(:,46))).*power(CONSTANTS(:,44), CONSTANTS(:,133))).*0.450000;
    ALGEBRAIC(:,107) =  (1.00000./CONSTANTS(:,110)).*log(CONSTANTS(:,54)./STATES(:,18));
    ALGEBRAIC(:,108) =  ( ( ( ( CONSTANTS(:,128).*CONSTANTS(:,116)).*power(STATES(:,34), 3.00000)).*STATES(:,32)).*STATES(:,33)).*(STATES(:,21) - ALGEBRAIC(:,107));
    ALGEBRAIC(:,110) =  ( CONSTANTS(:,128).*CONSTANTS(:,65)).*(STATES(:,21) - ALGEBRAIC(:,107));
    ALGEBRAIC(:,73) = (( ( ( CONSTANTS(:,128).*CONSTANTS(:,75)).*ALGEBRAIC(:,70)).*CONSTANTS(:,53))./(1.00000+power(CONSTANTS(:,118)./STATES(:,18), 4.00000)))./(CONSTANTS(:,53)+CONSTANTS(:,76));
    ALGEBRAIC(:,112) =  ( ( ( CONSTANTS(:,128).*CONSTANTS(:,120)).*power(STATES(:,36), 3.00000)).*STATES(:,35)).*(STATES(:,21) - ALGEBRAIC(:,107));
    ALGEBRAIC(:,116) = ((((ALGEBRAIC(:,108)+ALGEBRAIC(:,110))+ 3.00000.*ALGEBRAIC(:,69))+ 3.00000.*ALGEBRAIC(:,73))+ALGEBRAIC(:,47))+ALGEBRAIC(:,112);
    RATES(:,18) = (((  - ALGEBRAIC(:,116).*CONSTANTS(:,30))./( CONSTANTS(:,131).*CONSTANTS(:,32))+ (CONSTANTS(:,37)./CONSTANTS(:,131)).*(STATES(:,17) - STATES(:,18)))+ (CONSTANTS(:,38)./CONSTANTS(:,131)).*(STATES(:,29) - STATES(:,18))) - RATES(:,16);
    ALGEBRAIC(:,91) = ALGEBRAIC(:,83)+ALGEBRAIC(:,90);
    ALGEBRAIC(:,48) =  CONSTANTS(:,55).*(STATES(:,21) - CONSTANTS(:,126));
    ALGEBRAIC(:,50) =  (( CONSTANTS(:,43).*CONSTANTS(:,56))./(1.00000+CONSTANTS(:,57)./STATES(:,13))).*(STATES(:,21) - CONSTANTS(:,126));
    ALGEBRAIC(:,52) =  (( CONSTANTS(:,128).*CONSTANTS(:,56))./(1.00000+CONSTANTS(:,57)./STATES(:,14))).*(STATES(:,21) - CONSTANTS(:,126));
    ALGEBRAIC(:,53) = ALGEBRAIC(:,50)+ALGEBRAIC(:,52);
    ALGEBRAIC(:,79) = ALGEBRAIC(:,53)+ALGEBRAIC(:,48);
    ALGEBRAIC(:,117) = ALGEBRAIC(:,115)+ALGEBRAIC(:,116);
    ALGEBRAIC(:,42) =  (( ( ( ( CONSTANTS(:,107).*STATES(:,21)).*CONSTANTS(:,32)).*CONSTANTS(:,110)).*( ( 0.750000.*STATES(:,26)).*exp( STATES(:,21).*CONSTANTS(:,110)) -  0.750000.*CONSTANTS(:,53)))./(exp( STATES(:,21).*CONSTANTS(:,110)) - 1.00000)).*CONSTANTS(:,45);
    ALGEBRAIC(:,43) =  ( ( ( ( ALGEBRAIC(:,42).*STATES(:,22)).*STATES(:,23)).*( CONSTANTS(:,51).*(CONSTANTS(:,47)+(1.00000 - STATES(:,24)))+ CONSTANTS(:,109).*(CONSTANTS(:,46)+(1.00000 - STATES(:,25))))).*power(CONSTANTS(:,44), CONSTANTS(:,133))).*0.450000;
    ALGEBRAIC(:,92) =  (1.00000./CONSTANTS(:,110)).*log(CONSTANTS(:,53)./STATES(:,26));
    ALGEBRAIC(:,93) =  1.00000.*(1.02000./(1.00000+exp( 0.238500.*((STATES(:,21) - ALGEBRAIC(:,92)) - 59.2150))));
    ALGEBRAIC(:,94) =  1.00000.*(( 0.491240.*exp( 0.0803200.*((STATES(:,21)+5.47600) - ALGEBRAIC(:,92)))+exp( 0.0617500.*((STATES(:,21) - ALGEBRAIC(:,92)) - 594.310)))./(1.00000+exp(  - 0.514300.*((STATES(:,21) - ALGEBRAIC(:,92))+4.75300))));
    ALGEBRAIC(:,95) = ALGEBRAIC(:,93)./(ALGEBRAIC(:,93)+ALGEBRAIC(:,94));
    ALGEBRAIC(:,96) =  ( CONSTANTS(:,111).*ALGEBRAIC(:,95)).*(STATES(:,21) - ALGEBRAIC(:,92));
    ALGEBRAIC(:,54) = 1.00000./(1.00000+exp(7.48800 - STATES(:,21)./5.98000));
    ALGEBRAIC(:,97) =  ( ( CONSTANTS(:,43).*CONSTANTS(:,59)).*ALGEBRAIC(:,54)).*(STATES(:,21) - ALGEBRAIC(:,92));
    ALGEBRAIC(:,98) =  ( ( CONSTANTS(:,128).*CONSTANTS(:,59)).*ALGEBRAIC(:,54)).*(STATES(:,21) - ALGEBRAIC(:,92));
    ALGEBRAIC(:,99) = ALGEBRAIC(:,97)+ALGEBRAIC(:,98);
    ALGEBRAIC(:,55) = 1.00000./(1.00000+exp((STATES(:,21)+74.0000)./24.0000));
    ALGEBRAIC(:,100) =  ( ( CONSTANTS(:,112).*STATES(:,27)).*ALGEBRAIC(:,55)).*(STATES(:,21) - ALGEBRAIC(:,92));
    ALGEBRAIC(:,56) =  (1.00000./CONSTANTS(:,110)).*log((CONSTANTS(:,53)+ CONSTANTS(:,62).*CONSTANTS(:,54))./(STATES(:,26)+ CONSTANTS(:,62).*STATES(:,29)));
    ALGEBRAIC(:,57) =  ( ( CONSTANTS(:,43).*CONSTANTS(:,113)).*power(STATES(:,28), 2.00000)).*(STATES(:,21) - ALGEBRAIC(:,56));
    ALGEBRAIC(:,58) =  ( ( CONSTANTS(:,128).*CONSTANTS(:,114)).*power(STATES(:,28), 2.00000)).*(STATES(:,21) - ALGEBRAIC(:,56));
    ALGEBRAIC(:,59) = ALGEBRAIC(:,57)+ALGEBRAIC(:,58);
    ALGEBRAIC(:,101) =  ( ( CONSTANTS(:,115).*STATES(:,30)).*STATES(:,31)).*(STATES(:,21) - ALGEBRAIC(:,92));
    ALGEBRAIC(:,74) = ALGEBRAIC(:,72)+ALGEBRAIC(:,73);
    ALGEBRAIC(:,102) =  ( ( CONSTANTS(:,121).*STATES(:,37)).*STATES(:,38)).*(STATES(:,21) - ALGEBRAIC(:,92));
    ALGEBRAIC(:,114) = ((((((ALGEBRAIC(:,102)+ALGEBRAIC(:,100))+ALGEBRAIC(:,59))+ALGEBRAIC(:,96)) -  2.00000.*ALGEBRAIC(:,74))+ALGEBRAIC(:,43))+ALGEBRAIC(:,99))+ALGEBRAIC(:,101);
    ALGEBRAIC(:,118) = ((ALGEBRAIC(:,117)+ALGEBRAIC(:,79))+ALGEBRAIC(:,91))+ALGEBRAIC(:,114);
    ALGEBRAIC(:,119) = piecewise({(VOI - CONSTANTS(:,101)) -  CONSTANTS(:,102).*floor((VOI - CONSTANTS(:,101))./CONSTANTS(:,102))<CONSTANTS(:,100), 1.00000 }, 0.00000);
    ALGEBRAIC(:,120) =  ALGEBRAIC(:,119).*CONSTANTS(:,99);
    RATES(:,21) =  - (ALGEBRAIC(:,118)+ALGEBRAIC(:,120));
   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(:,10) = 1.00000./(1.00000+exp((STATES(:,21)+91.0000)./6.10000));
    ALGEBRAIC(:,1) = 1.00000./(1.00000+exp( - ((STATES(:,21)+ 3.00000.*CONSTANTS(:,24))+9.00000)./6.00000));
    ALGEBRAIC(:,15) = ( ALGEBRAIC(:,1).*(1.00000 - exp( - ((STATES(:,21)+ 3.00000.*CONSTANTS(:,24))+9.00000)./6.00000)))./( 0.0350000.*((STATES(:,21)+ 3.00000.*CONSTANTS(:,24))+9.00000));
    ALGEBRAIC(:,2) = 1.00000./(1.00000+exp(((STATES(:,21)+ 3.00000.*CONSTANTS(:,24))+30.0000)./7.00000))+0.200000./(1.00000+exp(((50.0000 - STATES(:,21)) -  3.00000.*CONSTANTS(:,24))./20.0000));
    ALGEBRAIC(:,16) = 1.00000./( 0.0197000.*exp( - power( 0.0337000.*((STATES(:,21)+ 3.00000.*CONSTANTS(:,24))+25.0000), 2.00000))+0.0200000);
    ALGEBRAIC(:,3) = 1.00000./(1.00000+exp( - (STATES(:,21)+10.0000)./5.00000));
    ALGEBRAIC(:,17) = ( ((550.000)./(1.00000+exp(( - 22.0000 - STATES(:,21))./9.00000))).*6.00000)./(1.00000+exp((STATES(:,21)+11.0000)./9.00000))+230.000./(1.00000+exp((STATES(:,21)+40.0000)./20.0000));
    ALGEBRAIC(:,4) = 1.00000./(1.00000+exp( - ((STATES(:,21)+ 40.0000.*CONSTANTS(:,24))+3.80000)./14.2500));
    ALGEBRAIC(:,18) = 990.100./(1.00000+exp( - ((STATES(:,21)+ 40.0000.*CONSTANTS(:,24))+2.43600)./14.1200));
    ALGEBRAIC(:,5) = 1.00000./(1.00000+exp((STATES(:,21)+6.00000)./ - 8.60000));
    ALGEBRAIC(:,19) = 9.00000./(1.00000+exp((STATES(:,21)+5.00000)./12.0000))+0.500000;
    ALGEBRAIC(:,6) = 1.00000./(1.00000+exp((STATES(:,21)+7.50000)./10.0000));
    ALGEBRAIC(:,20) = 590.000./(1.00000+exp((STATES(:,21)+60.0000)./10.0000))+3050.00;
    ALGEBRAIC(:,9) = 1.00000./power(1.00000+exp( - (56.8600+STATES(:,21))./9.03000), 2.00000);
    ALGEBRAIC(:,23) =  0.129200.*exp( - power((STATES(:,21)+45.7900)./15.5400, 2.00000))+ 0.0648700.*exp( - power((STATES(:,21) - 4.82300)./51.1200, 2.00000));
    ALGEBRAIC(:,11) = ( 0.320000.*(STATES(:,21)+47.1300))./(1.00000 - exp(  - 0.100000.*(STATES(:,21)+47.1300)));
    ALGEBRAIC(:,24) =  0.0800000.*exp( - STATES(:,21)./11.0000);
    ALGEBRAIC(:,12) = 1.00000./(1.00000+exp( - (STATES(:,21)+1.00000)./11.0000));
    ALGEBRAIC(:,25) =  3.50000.*exp( - power(STATES(:,21)./30.0000, 2.00000))+1.50000;
    ALGEBRAIC(:,13) = 1.00000./(1.00000+exp((STATES(:,21)+40.5000)./11.5000));
    ALGEBRAIC(:,26) =  25.6350.*exp( - power((STATES(:,21)+52.4500)./15.8827, 2.00000))+24.1400;
    ALGEBRAIC(:,28) = 1.00000./power(1.00000+exp((STATES(:,21)+71.5500)./7.43000), 2.00000);
    ALGEBRAIC(:,7) =  piecewise({STATES(:,21)>= - 40.0000, 0.00000 },  0.0570000.*exp( - (STATES(:,21)+80.0000)./6.80000)).*1.00000;
    ALGEBRAIC(:,21) =  piecewise({STATES(:,21)>= - 40.0000, 0.770000./( 0.130000.*(1.00000+exp( - (STATES(:,21)+10.6600)./11.1000))) },  2.70000.*exp( 0.0790000.*STATES(:,21))+ ( 3.10000.*power(10.0000, 5.00000)).*exp( 0.348500.*STATES(:,21))).*1.00000;
    ALGEBRAIC(:,31) = 1.00000./(ALGEBRAIC(:,7)+ALGEBRAIC(:,21));
    ALGEBRAIC(:,29) = 1.00000./power(1.00000+exp((STATES(:,21)+71.5500)./7.43000), 2.00000);
    ALGEBRAIC(:,8) =  piecewise({STATES(:,21)>= - 40.0000, 0.00000 },  (( ( (  - 2.54280.*power(10.0000, 4.00000)).*exp( 0.244400.*STATES(:,21)) -  ( 6.94800.*power(10.0000,  - 6.00000)).*exp(  - 0.0439100.*STATES(:,21))).*(STATES(:,21)+37.7800))./(1.00000+exp( 0.311000.*(STATES(:,21)+79.2300)))).*1.00000).*1.00000;
    ALGEBRAIC(:,22) =  piecewise({STATES(:,21)>= - 40.0000, ( 0.600000.*exp( 0.0570000.*STATES(:,21)))./(1.00000+exp(  - 0.100000.*(STATES(:,21)+32.0000))) }, ( 0.0242400.*exp(  - 0.0105200.*STATES(:,21)))./(1.00000+exp(  - 0.137800.*(STATES(:,21)+40.1400)))).*1.00000;
    ALGEBRAIC(:,32) = 1.00000./(ALGEBRAIC(:,8)+ALGEBRAIC(:,22));
    ALGEBRAIC(:,14) = ((1.00000 - STATES(:,41)) - STATES(:,40)) - STATES(:,39);
    ALGEBRAIC(:,27) = CONSTANTS(:,89) - (CONSTANTS(:,89) - CONSTANTS(:,90))./(1.00000+power(CONSTANTS(:,93)./STATES(:,19), 2.50000));
    ALGEBRAIC(:,30) =  CONSTANTS(:,95).*ALGEBRAIC(:,27);
    ALGEBRAIC(:,33) = CONSTANTS(:,123)./ALGEBRAIC(:,27);
    ALGEBRAIC(:,34) = RATES(:,6)+RATES(:,4);
    ALGEBRAIC(:,81) =  ((1.00000./CONSTANTS(:,110))./2.00000).*log(CONSTANTS(:,52)./STATES(:,13));
    ALGEBRAIC(:,82) =  ( CONSTANTS(:,43).*CONSTANTS(:,42)).*(STATES(:,21) - ALGEBRAIC(:,81));
    ALGEBRAIC(:,37) =  (( ( ( CONSTANTS(:,106).*4.00000).*( ( STATES(:,21).*CONSTANTS(:,32)).*CONSTANTS(:,110))).*( ( 0.341000.*STATES(:,13)).*exp( ( 2.00000.*STATES(:,21)).*CONSTANTS(:,110)) -  0.341000.*CONSTANTS(:,52)))./(exp( ( 2.00000.*STATES(:,21)).*CONSTANTS(:,110)) - 1.00000)).*CONSTANTS(:,45);
    ALGEBRAIC(:,38) =  ( ( ( ( ( CONSTANTS(:,51).*ALGEBRAIC(:,37)).*STATES(:,22)).*STATES(:,23)).*((1.00000 - STATES(:,24))+CONSTANTS(:,47))).*power(CONSTANTS(:,44), CONSTANTS(:,133))).*0.450000;
    ALGEBRAIC(:,60) = 1.00000./(1.00000+power(CONSTANTS(:,67)./STATES(:,13), 2.00000));
    ALGEBRAIC(:,62) =  ( exp( ( CONSTANTS(:,74).*STATES(:,21)).*CONSTANTS(:,110)).*power(STATES(:,17), 3.00000)).*CONSTANTS(:,52);
    ALGEBRAIC(:,63) =  ( exp( ( (CONSTANTS(:,74) - 1.00000).*STATES(:,21)).*CONSTANTS(:,110)).*power(CONSTANTS(:,54), 3.00000)).*STATES(:,13);
    ALGEBRAIC(:,64) = ((( ( CONSTANTS(:,68).*power(CONSTANTS(:,54), 3.00000)).*(1.00000+power(STATES(:,17)./CONSTANTS(:,70), 3.00000))+ ( power(CONSTANTS(:,71), 3.00000).*STATES(:,13)).*(1.00000+STATES(:,13)./CONSTANTS(:,68)))+ CONSTANTS(:,69).*power(STATES(:,17), 3.00000))+ power(STATES(:,17), 3.00000).*CONSTANTS(:,52))+ power(CONSTANTS(:,54), 3.00000).*STATES(:,13);
    ALGEBRAIC(:,65) = (( ( ( ( CONSTANTS(:,43).*CONSTANTS(:,117)).*power(CONSTANTS(:,72), CONSTANTS(:,133))).*ALGEBRAIC(:,60)).*(ALGEBRAIC(:,62) - ALGEBRAIC(:,63)))./ALGEBRAIC(:,64))./(1.00000+ CONSTANTS(:,73).*exp( ( (CONSTANTS(:,74) - 1.00000).*STATES(:,21)).*CONSTANTS(:,110)));
    ALGEBRAIC(:,75) = power( STATES(:,13).*1.00000, 1.60000);
    ALGEBRAIC(:,76) = ( ( ( CONSTANTS(:,43).*power(CONSTANTS(:,83), CONSTANTS(:,133))).*CONSTANTS(:,81)).*ALGEBRAIC(:,75))./(ALGEBRAIC(:,75)+CONSTANTS(:,105));
    ALGEBRAIC(:,83) = ((ALGEBRAIC(:,38)+ALGEBRAIC(:,82))+ALGEBRAIC(:,76)) -  2.00000.*ALGEBRAIC(:,65);
    ALGEBRAIC(:,84) =  ( CONSTANTS(:,98).*STATES(:,40)).*(STATES(:,19) - STATES(:,13));
    ALGEBRAIC(:,86) =  ( (1.00000+ 0.250000.*CONSTANTS(:,35)).*(STATES(:,19) - STATES(:,13))).*CONSTANTS(:,87);
    ALGEBRAIC(:,36) = (((((RATES(:,11)+RATES(:,9))+RATES(:,10))+RATES(:,1))+RATES(:,2))+RATES(:,3))+RATES(:,8);
    ALGEBRAIC(:,88) = ( ( power(CONSTANTS(:,91), CONSTANTS(:,133)).*CONSTANTS(:,92)).*(power(STATES(:,12)./CONSTANTS(:,122), CONSTANTS(:,94)) - power(STATES(:,19)./CONSTANTS(:,88), CONSTANTS(:,94))))./((1.00000+power(STATES(:,12)./CONSTANTS(:,122), CONSTANTS(:,94)))+power(STATES(:,19)./CONSTANTS(:,88), CONSTANTS(:,94)));
    ALGEBRAIC(:,35) = RATES(:,7)+RATES(:,5);
    ALGEBRAIC(:,85) =  ((1.00000./CONSTANTS(:,110))./2.00000).*log(CONSTANTS(:,52)./STATES(:,14));
    ALGEBRAIC(:,87) =  ( CONSTANTS(:,128).*CONSTANTS(:,42)).*(STATES(:,21) - ALGEBRAIC(:,85));
    ALGEBRAIC(:,39) =  (( ( ( CONSTANTS(:,106).*4.00000).*( ( STATES(:,21).*CONSTANTS(:,32)).*CONSTANTS(:,110))).*( ( 0.341000.*STATES(:,14)).*exp( ( 2.00000.*STATES(:,21)).*CONSTANTS(:,110)) -  0.341000.*CONSTANTS(:,52)))./(exp( ( 2.00000.*STATES(:,21)).*CONSTANTS(:,110)) - 1.00000)).*CONSTANTS(:,45);
    ALGEBRAIC(:,40) =  ( ( ( ( ( CONSTANTS(:,109).*ALGEBRAIC(:,39)).*STATES(:,22)).*STATES(:,23)).*((1.00000 - STATES(:,25))+CONSTANTS(:,46))).*power(CONSTANTS(:,44), CONSTANTS(:,133))).*0.450000;
    ALGEBRAIC(:,61) = 1.00000./(1.00000+power(CONSTANTS(:,67)./STATES(:,14), 2.00000));
    ALGEBRAIC(:,66) =  ( exp( ( CONSTANTS(:,74).*STATES(:,21)).*CONSTANTS(:,110)).*power(STATES(:,18), 3.00000)).*CONSTANTS(:,52);
    ALGEBRAIC(:,67) =  ( exp( ( (CONSTANTS(:,74) - 1.00000).*STATES(:,21)).*CONSTANTS(:,110)).*power(CONSTANTS(:,54), 3.00000)).*STATES(:,14);
    ALGEBRAIC(:,68) = ((( ( CONSTANTS(:,68).*power(CONSTANTS(:,54), 3.00000)).*(1.00000+power(STATES(:,18)./CONSTANTS(:,70), 3.00000))+ ( power(CONSTANTS(:,71), 3.00000).*STATES(:,14)).*(1.00000+STATES(:,14)./CONSTANTS(:,68)))+ CONSTANTS(:,69).*power(STATES(:,18), 3.00000))+ power(STATES(:,18), 3.00000).*CONSTANTS(:,52))+ power(CONSTANTS(:,54), 3.00000).*STATES(:,14);
    ALGEBRAIC(:,69) = (( ( ( ( CONSTANTS(:,128).*CONSTANTS(:,117)).*power(CONSTANTS(:,72), CONSTANTS(:,133))).*ALGEBRAIC(:,61)).*(ALGEBRAIC(:,66) - ALGEBRAIC(:,67)))./ALGEBRAIC(:,68))./(1.00000+ CONSTANTS(:,73).*exp( ( (CONSTANTS(:,74) - 1.00000).*STATES(:,21)).*CONSTANTS(:,110)));
    ALGEBRAIC(:,77) = power( STATES(:,14).*1.00000, 1.60000);
    ALGEBRAIC(:,78) = ( ( ( CONSTANTS(:,128).*power(CONSTANTS(:,83), CONSTANTS(:,133))).*CONSTANTS(:,81)).*ALGEBRAIC(:,77))./(ALGEBRAIC(:,77)+CONSTANTS(:,125));
    ALGEBRAIC(:,90) = ((ALGEBRAIC(:,40)+ALGEBRAIC(:,87))+ALGEBRAIC(:,78)) -  2.00000.*ALGEBRAIC(:,69);
    ALGEBRAIC(:,44) =  (( ( ( ( CONSTANTS(:,108).*STATES(:,21)).*CONSTANTS(:,32)).*CONSTANTS(:,110)).*( ( 0.750000.*STATES(:,17)).*exp( STATES(:,21).*CONSTANTS(:,110)) -  0.750000.*CONSTANTS(:,54)))./(exp( STATES(:,21).*CONSTANTS(:,110)) - 1.00000)).*CONSTANTS(:,45);
    ALGEBRAIC(:,45) =  ( ( ( ( ( CONSTANTS(:,51).*ALGEBRAIC(:,44)).*STATES(:,22)).*STATES(:,23)).*((1.00000 - STATES(:,24))+CONSTANTS(:,47))).*power(CONSTANTS(:,44), CONSTANTS(:,133))).*0.450000;
    ALGEBRAIC(:,103) =  (1.00000./CONSTANTS(:,110)).*log(CONSTANTS(:,54)./STATES(:,17));
    ALGEBRAIC(:,104) =  ( ( ( ( CONSTANTS(:,43).*CONSTANTS(:,116)).*power(STATES(:,34), 3.00000)).*STATES(:,32)).*STATES(:,33)).*(STATES(:,21) - ALGEBRAIC(:,103));
    ALGEBRAIC(:,105) =  ( CONSTANTS(:,43).*CONSTANTS(:,65)).*(STATES(:,21) - ALGEBRAIC(:,103));
    ALGEBRAIC(:,70) = 1.00000./((1.00000+ 0.124500.*exp( (  - 0.100000.*STATES(:,21)).*CONSTANTS(:,110)))+ ( 0.0365000.*CONSTANTS(:,119)).*exp(  - STATES(:,21).*CONSTANTS(:,110)));
    ALGEBRAIC(:,72) = (( ( ( CONSTANTS(:,43).*CONSTANTS(:,75)).*ALGEBRAIC(:,70)).*CONSTANTS(:,53))./(1.00000+power(CONSTANTS(:,118)./STATES(:,17), 4.00000)))./(CONSTANTS(:,53)+CONSTANTS(:,76));
    ALGEBRAIC(:,106) =  ( ( ( CONSTANTS(:,43).*CONSTANTS(:,120)).*power(STATES(:,36), 3.00000)).*STATES(:,35)).*(STATES(:,21) - ALGEBRAIC(:,103));
    ALGEBRAIC(:,115) = ((((ALGEBRAIC(:,104)+ALGEBRAIC(:,105))+ 3.00000.*ALGEBRAIC(:,65))+ 3.00000.*ALGEBRAIC(:,72))+ALGEBRAIC(:,45))+ALGEBRAIC(:,106);
    ALGEBRAIC(:,46) =  (( ( ( ( CONSTANTS(:,108).*STATES(:,21)).*CONSTANTS(:,32)).*CONSTANTS(:,110)).*( ( 0.750000.*STATES(:,18)).*exp( STATES(:,21).*CONSTANTS(:,110)) -  0.750000.*CONSTANTS(:,54)))./(exp( STATES(:,21).*CONSTANTS(:,110)) - 1.00000)).*CONSTANTS(:,45);
    ALGEBRAIC(:,47) =  ( ( ( ( ( CONSTANTS(:,109).*ALGEBRAIC(:,46)).*STATES(:,22)).*STATES(:,23)).*((1.00000 - STATES(:,25))+CONSTANTS(:,46))).*power(CONSTANTS(:,44), CONSTANTS(:,133))).*0.450000;
    ALGEBRAIC(:,107) =  (1.00000./CONSTANTS(:,110)).*log(CONSTANTS(:,54)./STATES(:,18));
    ALGEBRAIC(:,108) =  ( ( ( ( CONSTANTS(:,128).*CONSTANTS(:,116)).*power(STATES(:,34), 3.00000)).*STATES(:,32)).*STATES(:,33)).*(STATES(:,21) - ALGEBRAIC(:,107));
    ALGEBRAIC(:,110) =  ( CONSTANTS(:,128).*CONSTANTS(:,65)).*(STATES(:,21) - ALGEBRAIC(:,107));
    ALGEBRAIC(:,73) = (( ( ( CONSTANTS(:,128).*CONSTANTS(:,75)).*ALGEBRAIC(:,70)).*CONSTANTS(:,53))./(1.00000+power(CONSTANTS(:,118)./STATES(:,18), 4.00000)))./(CONSTANTS(:,53)+CONSTANTS(:,76));
    ALGEBRAIC(:,112) =  ( ( ( CONSTANTS(:,128).*CONSTANTS(:,120)).*power(STATES(:,36), 3.00000)).*STATES(:,35)).*(STATES(:,21) - ALGEBRAIC(:,107));
    ALGEBRAIC(:,116) = ((((ALGEBRAIC(:,108)+ALGEBRAIC(:,110))+ 3.00000.*ALGEBRAIC(:,69))+ 3.00000.*ALGEBRAIC(:,73))+ALGEBRAIC(:,47))+ALGEBRAIC(:,112);
    ALGEBRAIC(:,91) = ALGEBRAIC(:,83)+ALGEBRAIC(:,90);
    ALGEBRAIC(:,48) =  CONSTANTS(:,55).*(STATES(:,21) - CONSTANTS(:,126));
    ALGEBRAIC(:,50) =  (( CONSTANTS(:,43).*CONSTANTS(:,56))./(1.00000+CONSTANTS(:,57)./STATES(:,13))).*(STATES(:,21) - CONSTANTS(:,126));
    ALGEBRAIC(:,52) =  (( CONSTANTS(:,128).*CONSTANTS(:,56))./(1.00000+CONSTANTS(:,57)./STATES(:,14))).*(STATES(:,21) - CONSTANTS(:,126));
    ALGEBRAIC(:,53) = ALGEBRAIC(:,50)+ALGEBRAIC(:,52);
    ALGEBRAIC(:,79) = ALGEBRAIC(:,53)+ALGEBRAIC(:,48);
    ALGEBRAIC(:,117) = ALGEBRAIC(:,115)+ALGEBRAIC(:,116);
    ALGEBRAIC(:,42) =  (( ( ( ( CONSTANTS(:,107).*STATES(:,21)).*CONSTANTS(:,32)).*CONSTANTS(:,110)).*( ( 0.750000.*STATES(:,26)).*exp( STATES(:,21).*CONSTANTS(:,110)) -  0.750000.*CONSTANTS(:,53)))./(exp( STATES(:,21).*CONSTANTS(:,110)) - 1.00000)).*CONSTANTS(:,45);
    ALGEBRAIC(:,43) =  ( ( ( ( ALGEBRAIC(:,42).*STATES(:,22)).*STATES(:,23)).*( CONSTANTS(:,51).*(CONSTANTS(:,47)+(1.00000 - STATES(:,24)))+ CONSTANTS(:,109).*(CONSTANTS(:,46)+(1.00000 - STATES(:,25))))).*power(CONSTANTS(:,44), CONSTANTS(:,133))).*0.450000;
    ALGEBRAIC(:,92) =  (1.00000./CONSTANTS(:,110)).*log(CONSTANTS(:,53)./STATES(:,26));
    ALGEBRAIC(:,93) =  1.00000.*(1.02000./(1.00000+exp( 0.238500.*((STATES(:,21) - ALGEBRAIC(:,92)) - 59.2150))));
    ALGEBRAIC(:,94) =  1.00000.*(( 0.491240.*exp( 0.0803200.*((STATES(:,21)+5.47600) - ALGEBRAIC(:,92)))+exp( 0.0617500.*((STATES(:,21) - ALGEBRAIC(:,92)) - 594.310)))./(1.00000+exp(  - 0.514300.*((STATES(:,21) - ALGEBRAIC(:,92))+4.75300))));
    ALGEBRAIC(:,95) = ALGEBRAIC(:,93)./(ALGEBRAIC(:,93)+ALGEBRAIC(:,94));
    ALGEBRAIC(:,96) =  ( CONSTANTS(:,111).*ALGEBRAIC(:,95)).*(STATES(:,21) - ALGEBRAIC(:,92));
    ALGEBRAIC(:,54) = 1.00000./(1.00000+exp(7.48800 - STATES(:,21)./5.98000));
    ALGEBRAIC(:,97) =  ( ( CONSTANTS(:,43).*CONSTANTS(:,59)).*ALGEBRAIC(:,54)).*(STATES(:,21) - ALGEBRAIC(:,92));
    ALGEBRAIC(:,98) =  ( ( CONSTANTS(:,128).*CONSTANTS(:,59)).*ALGEBRAIC(:,54)).*(STATES(:,21) - ALGEBRAIC(:,92));
    ALGEBRAIC(:,99) = ALGEBRAIC(:,97)+ALGEBRAIC(:,98);
    ALGEBRAIC(:,55) = 1.00000./(1.00000+exp((STATES(:,21)+74.0000)./24.0000));
    ALGEBRAIC(:,100) =  ( ( CONSTANTS(:,112).*STATES(:,27)).*ALGEBRAIC(:,55)).*(STATES(:,21) - ALGEBRAIC(:,92));
    ALGEBRAIC(:,56) =  (1.00000./CONSTANTS(:,110)).*log((CONSTANTS(:,53)+ CONSTANTS(:,62).*CONSTANTS(:,54))./(STATES(:,26)+ CONSTANTS(:,62).*STATES(:,29)));
    ALGEBRAIC(:,57) =  ( ( CONSTANTS(:,43).*CONSTANTS(:,113)).*power(STATES(:,28), 2.00000)).*(STATES(:,21) - ALGEBRAIC(:,56));
    ALGEBRAIC(:,58) =  ( ( CONSTANTS(:,128).*CONSTANTS(:,114)).*power(STATES(:,28), 2.00000)).*(STATES(:,21) - ALGEBRAIC(:,56));
    ALGEBRAIC(:,59) = ALGEBRAIC(:,57)+ALGEBRAIC(:,58);
    ALGEBRAIC(:,101) =  ( ( CONSTANTS(:,115).*STATES(:,30)).*STATES(:,31)).*(STATES(:,21) - ALGEBRAIC(:,92));
    ALGEBRAIC(:,74) = ALGEBRAIC(:,72)+ALGEBRAIC(:,73);
    ALGEBRAIC(:,102) =  ( ( CONSTANTS(:,121).*STATES(:,37)).*STATES(:,38)).*(STATES(:,21) - ALGEBRAIC(:,92));
    ALGEBRAIC(:,114) = ((((((ALGEBRAIC(:,102)+ALGEBRAIC(:,100))+ALGEBRAIC(:,59))+ALGEBRAIC(:,96)) -  2.00000.*ALGEBRAIC(:,74))+ALGEBRAIC(:,43))+ALGEBRAIC(:,99))+ALGEBRAIC(:,101);
    ALGEBRAIC(:,118) = ((ALGEBRAIC(:,117)+ALGEBRAIC(:,79))+ALGEBRAIC(:,91))+ALGEBRAIC(:,114);
    ALGEBRAIC(:,119) = piecewise({(VOI - CONSTANTS(:,101)) -  CONSTANTS(:,102).*floor((VOI - CONSTANTS(:,101))./CONSTANTS(:,102))<CONSTANTS(:,100), 1.00000 }, 0.00000);
    ALGEBRAIC(:,120) =  ALGEBRAIC(:,119).*CONSTANTS(:,99);
    ALGEBRAIC(:,41) = ALGEBRAIC(:,38)+ALGEBRAIC(:,40);
    ALGEBRAIC(:,49) = ALGEBRAIC(:,45)+ALGEBRAIC(:,47);
    ALGEBRAIC(:,51) = (ALGEBRAIC(:,41)+ALGEBRAIC(:,43))+ALGEBRAIC(:,49);
    ALGEBRAIC(:,71) = ALGEBRAIC(:,65)+ALGEBRAIC(:,69);
    ALGEBRAIC(:,80) = ALGEBRAIC(:,76)+ALGEBRAIC(:,78);
    ALGEBRAIC(:,89) = ALGEBRAIC(:,82)+ALGEBRAIC(:,87);
    ALGEBRAIC(:,109) = ALGEBRAIC(:,104)+ALGEBRAIC(:,108);
    ALGEBRAIC(:,111) = ALGEBRAIC(:,105)+ALGEBRAIC(:,110);
    ALGEBRAIC(:,113) = ALGEBRAIC(:,106)+ALGEBRAIC(:,112);
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