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 =0;
end
% There are a total of 33 entries in each of the rate and state variable arrays.
% There are a total of 80 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 (s)');
    LEGEND_CONSTANTS(:,1) = strpad('kf1 in component PI3K (second_order_rate_constant)');
    LEGEND_CONSTANTS(:,2) = strpad('kb1 in component PI3K (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,3) = strpad('kf2 in component PI3K (second_order_rate_constant)');
    LEGEND_CONSTANTS(:,4) = strpad('kb2 in component PI3K (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,5) = strpad('kf3 in component PI3K (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,6) = strpad('kb3 in component PI3K (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,7) = strpad('kf34 in component PI3K (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,8) = strpad('kb34 in component PI3K (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,9) = strpad('V4 in component PI3K (flux)');
    LEGEND_CONSTANTS(:,10) = strpad('k4 in component PI3K (nm)');
    LEGEND_CONSTANTS(:,11) = strpad('kf5 in component PI3K (second_order_rate_constant)');
    LEGEND_CONSTANTS(:,12) = strpad('kb5 in component PI3K (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,13) = strpad('kf6 in component PI3K (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,14) = strpad('kb6 in component PI3K (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,15) = strpad('kf7 in component PI3K (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,16) = strpad('kb7 in component PI3K (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,17) = strpad('kf8 in component PI3K (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,18) = strpad('kb8 in component PI3K (second_order_rate_constant)');
    LEGEND_CONSTANTS(:,19) = strpad('kf9 in component PI3K (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,20) = strpad('kb9 in component PI3K (second_order_rate_constant)');
    LEGEND_CONSTANTS(:,21) = strpad('V10 in component PI3K (flux)');
    LEGEND_CONSTANTS(:,22) = strpad('k10 in component PI3K (nm)');
    LEGEND_CONSTANTS(:,23) = strpad('kf23 in component PI3K (second_order_rate_constant)');
    LEGEND_CONSTANTS(:,24) = strpad('kb23 in component PI3K (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,25) = strpad('kf24 in component PI3K (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,26) = strpad('kb24 in component PI3K (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,27) = strpad('kf25 in component PI3K (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,28) = strpad('kb25 in component PI3K (second_order_rate_constant)');
    LEGEND_CONSTANTS(:,29) = strpad('V26 in component PI3K (flux)');
    LEGEND_CONSTANTS(:,30) = strpad('k26 in component PI3K (nm)');
    LEGEND_STATES(:,1) = strpad('R in component PI3K (nm)');
    LEGEND_STATES(:,2) = strpad('Shc in component PI3K (nm)');
    LEGEND_STATES(:,3) = strpad('PI3K in component PI3K (nm)');
    LEGEND_STATES(:,4) = strpad('HRG in component PI3K (nm)');
    LEGEND_STATES(:,5) = strpad('R_HRG in component PI3K (nm)');
    LEGEND_STATES(:,6) = strpad('R_HRG2 in component PI3K (nm)');
    LEGEND_STATES(:,7) = strpad('Internalisation in component PI3K (nm)');
    LEGEND_STATES(:,8) = strpad('RP in component PI3K (nm)');
    LEGEND_STATES(:,9) = strpad('R_Shc in component PI3K (nm)');
    LEGEND_STATES(:,10) = strpad('R_ShP in component PI3K (nm)');
    LEGEND_STATES(:,11) = strpad('ShP in component PI3K (nm)');
    LEGEND_STATES(:,12) = strpad('R_ShGS in component PI3K (nm)');
    LEGEND_STATES(:,13) = strpad('ShGS in component PI3K (nm)');
    LEGEND_STATES(:,14) = strpad('GS in component PI3K (nm)');
    LEGEND_STATES(:,15) = strpad('R_PI3K in component PI3K (nm)');
    LEGEND_STATES(:,16) = strpad('R_PI3Kstar in component PI3K (nm)');
    LEGEND_STATES(:,17) = strpad('PI3Kstar in component PI3K (nm)');
    LEGEND_CONSTANTS(:,31) = strpad('two in component PI3K (dimensionless)');
    LEGEND_STATES(:,18) = strpad('RasGTP in component RasGDPtoRasGTP (nm)');
    LEGEND_CONSTANTS(:,32) = strpad('kf11 in component RasGDPtoRasGTP (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,33) = strpad('k11 in component RasGDPtoRasGTP (nm)');
    LEGEND_CONSTANTS(:,34) = strpad('V12 in component RasGDPtoRasGTP (flux)');
    LEGEND_CONSTANTS(:,35) = strpad('k12 in component RasGDPtoRasGTP (nm)');
    LEGEND_STATES(:,19) = strpad('RasGDP in component RasGDPtoRasGTP (nm)');
    LEGEND_STATES(:,20) = strpad('Akt_PIPP in component Akt (nm)');
    LEGEND_STATES(:,21) = strpad('RAF_star in component RAF (nm)');
    LEGEND_CONSTANTS(:,36) = strpad('kf13 in component RAF (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,37) = strpad('k13 in component RAF (nm)');
    LEGEND_CONSTANTS(:,38) = strpad('kf14 in component RAF (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,39) = strpad('k14 in component RAF (nm)');
    LEGEND_CONSTANTS(:,40) = strpad('E in component RAF (nm)');
    LEGEND_STATES(:,22) = strpad('RAF in component RAF (nm)');
    LEGEND_STATES(:,23) = strpad('MEKP in component MEK (nm)');
    LEGEND_STATES(:,24) = strpad('MEKPP in component MEK (nm)');
    LEGEND_CONSTANTS(:,41) = strpad('kf27 in component Akt (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,42) = strpad('k27 in component Akt (nm)');
    LEGEND_CONSTANTS(:,43) = strpad('V28 in component Akt (flux)');
    LEGEND_CONSTANTS(:,44) = strpad('k28 in component Akt (nm)');
    LEGEND_CONSTANTS(:,45) = strpad('kf29 in component Akt (second_order_rate_constant)');
    LEGEND_CONSTANTS(:,46) = strpad('kb29 in component Akt (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,47) = strpad('V30 in component Akt (flux)');
    LEGEND_CONSTANTS(:,48) = strpad('k30 in component Akt (nm)');
    LEGEND_CONSTANTS(:,49) = strpad('kf31 in component Akt (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,50) = strpad('k31 in component Akt (nm)');
    LEGEND_CONSTANTS(:,51) = strpad('V32 in component Akt (flux)');
    LEGEND_CONSTANTS(:,52) = strpad('k32 in component Akt (nm)');
    LEGEND_CONSTANTS(:,53) = strpad('kf33 in component Akt (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,54) = strpad('k33 in component Akt (nm)');
    LEGEND_CONSTANTS(:,55) = strpad('k16 in component Akt (nm)');
    LEGEND_CONSTANTS(:,56) = strpad('k18 in component Akt (nm)');
    LEGEND_STATES(:,25) = strpad('P in component Akt (nm)');
    LEGEND_STATES(:,26) = strpad('PIP3 in component Akt (nm)');
    LEGEND_STATES(:,27) = strpad('Akt in component Akt (nm)');
    LEGEND_STATES(:,28) = strpad('Akt_PIP3 in component Akt (nm)');
    LEGEND_STATES(:,29) = strpad('Akt_PIP in component Akt (nm)');
    LEGEND_CONSTANTS(:,57) = strpad('PP2A in component Akt (nm)');
    LEGEND_CONSTANTS(:,58) = strpad('one in component Akt (dimensionless)');
    LEGEND_CONSTANTS(:,59) = strpad('PP2A in component MEK (nm)');
    LEGEND_STATES(:,30) = strpad('MEK in component MEK (nm)');
    LEGEND_CONSTANTS(:,60) = strpad('kf15 in component MEK (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,61) = strpad('k15 in component MEK (nm)');
    LEGEND_CONSTANTS(:,62) = strpad('kf16 in component MEK (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,63) = strpad('k16 in component MEK (nm)');
    LEGEND_CONSTANTS(:,64) = strpad('kf17 in component MEK (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,65) = strpad('k17 in component MEK (nm)');
    LEGEND_CONSTANTS(:,66) = strpad('kf18 in component MEK (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,67) = strpad('k18 in component MEK (nm)');
    LEGEND_CONSTANTS(:,68) = strpad('k31 in component MEK (nm)');
    LEGEND_CONSTANTS(:,69) = strpad('k33 in component MEK (nm)');
    LEGEND_CONSTANTS(:,70) = strpad('one in component MEK (dimensionless)');
    LEGEND_CONSTANTS(:,71) = strpad('MKP3 in component ERK (nm)');
    LEGEND_STATES(:,31) = strpad('ERK in component ERK (nm)');
    LEGEND_STATES(:,32) = strpad('ERKP in component ERK (nm)');
    LEGEND_STATES(:,33) = strpad('ERKPP in component ERK (nm)');
    LEGEND_CONSTANTS(:,72) = strpad('kf19 in component ERK (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,73) = strpad('k19 in component ERK (nm)');
    LEGEND_CONSTANTS(:,74) = strpad('kf20 in component ERK (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,75) = strpad('k20 in component ERK (nm)');
    LEGEND_CONSTANTS(:,76) = strpad('kf21 in component ERK (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,77) = strpad('k21 in component ERK (nm)');
    LEGEND_CONSTANTS(:,78) = strpad('kf22 in component ERK (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,79) = strpad('k22 in component ERK (nm)');
    LEGEND_CONSTANTS(:,80) = strpad('one in component ERK (dimensionless)');
    LEGEND_RATES(:,1) = strpad('d/dt R in component PI3K (nm)');
    LEGEND_RATES(:,4) = strpad('d/dt HRG in component PI3K (nm)');
    LEGEND_RATES(:,5) = strpad('d/dt R_HRG in component PI3K (nm)');
    LEGEND_RATES(:,6) = strpad('d/dt R_HRG2 in component PI3K (nm)');
    LEGEND_RATES(:,8) = strpad('d/dt RP in component PI3K (nm)');
    LEGEND_RATES(:,7) = strpad('d/dt Internalisation in component PI3K (nm)');
    LEGEND_RATES(:,9) = strpad('d/dt R_Shc in component PI3K (nm)');
    LEGEND_RATES(:,2) = strpad('d/dt Shc in component PI3K (nm)');
    LEGEND_RATES(:,10) = strpad('d/dt R_ShP in component PI3K (nm)');
    LEGEND_RATES(:,14) = strpad('d/dt GS in component PI3K (nm)');
    LEGEND_RATES(:,11) = strpad('d/dt ShP in component PI3K (nm)');
    LEGEND_RATES(:,12) = strpad('d/dt R_ShGS in component PI3K (nm)');
    LEGEND_RATES(:,13) = strpad('d/dt ShGS in component PI3K (nm)');
    LEGEND_RATES(:,15) = strpad('d/dt R_PI3K in component PI3K (nm)');
    LEGEND_RATES(:,3) = strpad('d/dt PI3K in component PI3K (nm)');
    LEGEND_RATES(:,16) = strpad('d/dt R_PI3Kstar in component PI3K (nm)');
    LEGEND_RATES(:,17) = strpad('d/dt PI3Kstar in component PI3K (nm)');
    LEGEND_RATES(:,19) = strpad('d/dt RasGDP in component RasGDPtoRasGTP (nm)');
    LEGEND_RATES(:,18) = strpad('d/dt RasGTP in component RasGDPtoRasGTP (nm)');
    LEGEND_RATES(:,22) = strpad('d/dt RAF in component RAF (nm)');
    LEGEND_RATES(:,21) = strpad('d/dt RAF_star in component RAF (nm)');
    LEGEND_RATES(:,25) = strpad('d/dt P in component Akt (nm)');
    LEGEND_RATES(:,26) = strpad('d/dt PIP3 in component Akt (nm)');
    LEGEND_RATES(:,27) = strpad('d/dt Akt in component Akt (nm)');
    LEGEND_RATES(:,28) = strpad('d/dt Akt_PIP3 in component Akt (nm)');
    LEGEND_RATES(:,29) = strpad('d/dt Akt_PIP in component Akt (nm)');
    LEGEND_RATES(:,20) = strpad('d/dt Akt_PIPP in component Akt (nm)');
    LEGEND_RATES(:,30) = strpad('d/dt MEK in component MEK (nm)');
    LEGEND_RATES(:,23) = strpad('d/dt MEKP in component MEK (nm)');
    LEGEND_RATES(:,24) = strpad('d/dt MEKPP in component MEK (nm)');
    LEGEND_RATES(:,31) = strpad('d/dt ERK in component ERK (nm)');
    LEGEND_RATES(:,33) = strpad('d/dt ERKPP in component ERK (nm)');
    LEGEND_RATES(:,32) = strpad('d/dt ERKP in component ERK (nm)');
    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.0012;
    CONSTANTS(:,2) = 0.00076;
    CONSTANTS(:,3) = 0.01;
    CONSTANTS(:,4) = 0.1;
    CONSTANTS(:,5) = 1;
    CONSTANTS(:,6) = 0.01;
    CONSTANTS(:,7) = 0.001;
    CONSTANTS(:,8) = 0;
    CONSTANTS(:,9) = 62.5;
    CONSTANTS(:,10) = 50;
    CONSTANTS(:,11) = 0.1;
    CONSTANTS(:,12) = 1;
    CONSTANTS(:,13) = 20;
    CONSTANTS(:,14) = 5;
    CONSTANTS(:,15) = 60;
    CONSTANTS(:,16) = 546;
    CONSTANTS(:,17) = 2040;
    CONSTANTS(:,18) = 15700;
    CONSTANTS(:,19) = 40.8;
    CONSTANTS(:,20) = 0;
    CONSTANTS(:,21) = 0.0154;
    CONSTANTS(:,22) = 340;
    CONSTANTS(:,23) = 0.1;
    CONSTANTS(:,24) = 2;
    CONSTANTS(:,25) = 9.85;
    CONSTANTS(:,26) = 0.0985;
    CONSTANTS(:,27) = 45.8;
    CONSTANTS(:,28) = 0.047;
    CONSTANTS(:,29) = 2620;
    CONSTANTS(:,30) = 3680;
    STATES(:,1) = 80;
    STATES(:,2) = 1000;
    STATES(:,3) = 10;
    STATES(:,4) = 10;
    STATES(:,5) = 0;
    STATES(:,6) = 0;
    STATES(:,7) = 0;
    STATES(:,8) = 0;
    STATES(:,9) = 0;
    STATES(:,10) = 0;
    STATES(:,11) = 0;
    STATES(:,12) = 0;
    STATES(:,13) = 0;
    STATES(:,14) = 10;
    STATES(:,15) = 0;
    STATES(:,16) = 0;
    STATES(:,17) = 0;
    CONSTANTS(:,31) = 2;
    STATES(:,18) = 0;
    CONSTANTS(:,32) = 0.222;
    CONSTANTS(:,33) = 0.181;
    CONSTANTS(:,34) = 0.289;
    CONSTANTS(:,35) = 0.0571;
    STATES(:,19) = 120;
    STATES(:,20) = 0.0;
    STATES(:,21) = 100;
    CONSTANTS(:,36) = 1.53;
    CONSTANTS(:,37) = 11.7;
    CONSTANTS(:,38) = 0.00673;
    CONSTANTS(:,39) = 8.07;
    CONSTANTS(:,40) = 7;
    STATES(:,22) = 0;
    STATES(:,23) = 0;
    STATES(:,24) = 0;
    CONSTANTS(:,41) = 16.9;
    CONSTANTS(:,42) = 39.1;
    CONSTANTS(:,43) = 17000;
    CONSTANTS(:,44) = 9.02;
    CONSTANTS(:,45) = 507;
    CONSTANTS(:,46) = 234;
    CONSTANTS(:,47) = 20000;
    CONSTANTS(:,48) = 80000;
    CONSTANTS(:,49) = 0.107;
    CONSTANTS(:,50) = 4.35;
    CONSTANTS(:,51) = 20000;
    CONSTANTS(:,52) = 80000;
    CONSTANTS(:,53) = 0.211;
    CONSTANTS(:,54) = 12;
    CONSTANTS(:,55) = 2200;
    CONSTANTS(:,56) = 60;
    STATES(:,25) = 800;
    STATES(:,26) = 0;
    STATES(:,27) = 10;
    STATES(:,28) = 0;
    STATES(:,29) = 0;
    CONSTANTS(:,57) = 11.4;
    CONSTANTS(:,58) = 1;
    CONSTANTS(:,59) = 11.4;
    STATES(:,30) = 120;
    CONSTANTS(:,60) = 3.5;
    CONSTANTS(:,61) = 317;
    CONSTANTS(:,62) = 0.058;
    CONSTANTS(:,63) = 2200;
    CONSTANTS(:,64) = 2.9;
    CONSTANTS(:,65) = 317;
    CONSTANTS(:,66) = 0.058;
    CONSTANTS(:,67) = 60;
    CONSTANTS(:,68) = 4.35;
    CONSTANTS(:,69) = 12;
    CONSTANTS(:,70) = 1;
    CONSTANTS(:,71) = 2.4;
    STATES(:,31) = 1000;
    STATES(:,32) = 0;
    STATES(:,33) = 0;
    CONSTANTS(:,72) = 9.5;
    CONSTANTS(:,73) = 146000;
    CONSTANTS(:,74) = 0.3;
    CONSTANTS(:,75) = 160;
    CONSTANTS(:,76) = 16;
    CONSTANTS(:,77) = 146000;
    CONSTANTS(:,78) = 0.27;
    CONSTANTS(:,79) = 60;
    CONSTANTS(:,80) = 1;
    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(:,1) =   - CONSTANTS(:,1).*STATES(:,1).*STATES(:,4)+ CONSTANTS(:,2).*STATES(:,5);
    RATES(:,4) =   - CONSTANTS(:,1).*STATES(:,1).*STATES(:,4)+ CONSTANTS(:,2).*STATES(:,5);
    RATES(:,5) = ( CONSTANTS(:,1).*STATES(:,1).*STATES(:,4) -  CONSTANTS(:,2).*STATES(:,5)) -  CONSTANTS(:,31).*( CONSTANTS(:,3).*STATES(:,5).*STATES(:,5) -  CONSTANTS(:,4).*STATES(:,6));
    RATES(:,6) = (( CONSTANTS(:,3).*STATES(:,5).*STATES(:,5) -  CONSTANTS(:,4).*STATES(:,6)) - ( CONSTANTS(:,5).*STATES(:,6) -  CONSTANTS(:,6).*STATES(:,8)))+( CONSTANTS(:,9).*STATES(:,8))./(CONSTANTS(:,10)+STATES(:,8));
    RATES(:,8) = (((((( CONSTANTS(:,5).*STATES(:,6) -  CONSTANTS(:,6).*STATES(:,8)) - ( CONSTANTS(:,9).*STATES(:,8))./(CONSTANTS(:,10)+STATES(:,8))) - ( CONSTANTS(:,11).*STATES(:,8).*STATES(:,2) -  CONSTANTS(:,12).*STATES(:,9)))+( CONSTANTS(:,17).*STATES(:,10) -  CONSTANTS(:,18).*STATES(:,13).*STATES(:,8))) - ( CONSTANTS(:,23).*STATES(:,8).*STATES(:,3) -  CONSTANTS(:,24).*STATES(:,15)))+( CONSTANTS(:,27).*STATES(:,16) -  CONSTANTS(:,28).*STATES(:,8).*STATES(:,17))) - ( CONSTANTS(:,7).*STATES(:,8) -  CONSTANTS(:,8).*STATES(:,7));
    RATES(:,7) =  CONSTANTS(:,7).*STATES(:,8) -  CONSTANTS(:,8).*STATES(:,7);
    RATES(:,9) = ( CONSTANTS(:,11).*STATES(:,8).*STATES(:,2) -  CONSTANTS(:,12).*STATES(:,9)) - ( CONSTANTS(:,13).*STATES(:,9) -  CONSTANTS(:,14).*STATES(:,10));
    RATES(:,2) =  - ( CONSTANTS(:,11).*STATES(:,8).*STATES(:,2) -  CONSTANTS(:,12).*STATES(:,9))+( CONSTANTS(:,21).*STATES(:,11))./(CONSTANTS(:,22)+STATES(:,11));
    RATES(:,10) = ( CONSTANTS(:,13).*STATES(:,9) -  CONSTANTS(:,14).*STATES(:,10)) - ( CONSTANTS(:,15).*STATES(:,10) -  CONSTANTS(:,16).*STATES(:,12));
    RATES(:,14) =  - ( CONSTANTS(:,15).*STATES(:,10) -  CONSTANTS(:,16).*STATES(:,12))+( CONSTANTS(:,19).*STATES(:,13) -  CONSTANTS(:,20).*STATES(:,14).*STATES(:,11));
    RATES(:,11) = ( CONSTANTS(:,19).*STATES(:,13) -  CONSTANTS(:,20).*STATES(:,14).*STATES(:,11)) - ( CONSTANTS(:,21).*STATES(:,11))./(CONSTANTS(:,22)+STATES(:,11));
    RATES(:,12) = ( CONSTANTS(:,15).*STATES(:,10) -  CONSTANTS(:,16).*STATES(:,12)) - ( CONSTANTS(:,17).*STATES(:,10) -  CONSTANTS(:,18).*STATES(:,13).*STATES(:,8));
    RATES(:,13) = ( CONSTANTS(:,17).*STATES(:,10) -  CONSTANTS(:,18).*STATES(:,13).*STATES(:,8)) - ( CONSTANTS(:,19).*STATES(:,13) -  CONSTANTS(:,20).*STATES(:,14).*STATES(:,11));
    RATES(:,15) = ( CONSTANTS(:,23).*STATES(:,8).*STATES(:,3) -  CONSTANTS(:,24).*STATES(:,15)) - ( CONSTANTS(:,25).*STATES(:,15) -  CONSTANTS(:,26).*STATES(:,16));
    RATES(:,3) =  - ( CONSTANTS(:,23).*STATES(:,8).*STATES(:,3) -  CONSTANTS(:,24).*STATES(:,15))+( CONSTANTS(:,29).*STATES(:,17))./(CONSTANTS(:,30)+STATES(:,17));
    RATES(:,16) = ( CONSTANTS(:,25).*STATES(:,15) -  CONSTANTS(:,26).*STATES(:,16)) - ( CONSTANTS(:,27).*STATES(:,16) -  CONSTANTS(:,28).*STATES(:,8).*STATES(:,17));
    RATES(:,17) = ( CONSTANTS(:,27).*STATES(:,16) -  CONSTANTS(:,28).*STATES(:,8).*STATES(:,17)) - ( CONSTANTS(:,29).*STATES(:,17))./(CONSTANTS(:,30)+STATES(:,17));
    RATES(:,19) =  - (( CONSTANTS(:,32).*STATES(:,13).*STATES(:,19))./(CONSTANTS(:,33)+STATES(:,19)))+( CONSTANTS(:,34).*STATES(:,18))./(CONSTANTS(:,35)+STATES(:,18));
    RATES(:,18) = ( CONSTANTS(:,32).*STATES(:,13).*STATES(:,19))./(CONSTANTS(:,33)+STATES(:,19)) - ( CONSTANTS(:,34).*STATES(:,18))./(CONSTANTS(:,35)+STATES(:,18));
    RATES(:,22) = ( CONSTANTS(:,38).*(STATES(:,20)+CONSTANTS(:,40)).*STATES(:,21))./(CONSTANTS(:,39)+STATES(:,21)) - ( CONSTANTS(:,36).*STATES(:,18).*STATES(:,22))./(CONSTANTS(:,37)+STATES(:,22));
    RATES(:,21) =  - (( CONSTANTS(:,38).*(STATES(:,20)+CONSTANTS(:,40)).*STATES(:,21))./(CONSTANTS(:,39)+STATES(:,21)))+( CONSTANTS(:,36).*STATES(:,18).*STATES(:,22))./(CONSTANTS(:,37)+STATES(:,22));
    RATES(:,25) = ( CONSTANTS(:,43).*STATES(:,26))./(CONSTANTS(:,44)+STATES(:,26)) - ( CONSTANTS(:,41).*STATES(:,17).*STATES(:,25))./(CONSTANTS(:,42)+STATES(:,25));
    RATES(:,26) = ( - (( CONSTANTS(:,43).*STATES(:,26))./(CONSTANTS(:,44)+STATES(:,26)))+( CONSTANTS(:,41).*STATES(:,17).*STATES(:,25))./(CONSTANTS(:,42)+STATES(:,25))) - ( CONSTANTS(:,45).*STATES(:,26).*STATES(:,27) -  CONSTANTS(:,46).*STATES(:,28));
    RATES(:,27) =  - ( CONSTANTS(:,45).*STATES(:,26).*STATES(:,27) -  CONSTANTS(:,46).*STATES(:,28));
    RATES(:,28) = (( CONSTANTS(:,45).*STATES(:,26).*STATES(:,27) -  CONSTANTS(:,46).*STATES(:,28)) - ( CONSTANTS(:,47).*STATES(:,28))./( CONSTANTS(:,48).*(CONSTANTS(:,58)+STATES(:,29)./CONSTANTS(:,52))+STATES(:,28)))+( CONSTANTS(:,49).*CONSTANTS(:,57).*STATES(:,29))./( CONSTANTS(:,50).*(CONSTANTS(:,58)+STATES(:,23)./CONSTANTS(:,55)+STATES(:,24)./CONSTANTS(:,56)+STATES(:,20)./CONSTANTS(:,54))+STATES(:,29));
    RATES(:,29) = ((( CONSTANTS(:,47).*STATES(:,28))./( CONSTANTS(:,48).*(CONSTANTS(:,58)+STATES(:,29)./CONSTANTS(:,52))+STATES(:,28)) - ( CONSTANTS(:,49).*CONSTANTS(:,57).*STATES(:,29))./( CONSTANTS(:,50).*(CONSTANTS(:,58)+STATES(:,23)./CONSTANTS(:,55)+STATES(:,24)./CONSTANTS(:,56)+STATES(:,20)./CONSTANTS(:,54))+STATES(:,29))) - ( CONSTANTS(:,51).*STATES(:,29))./( CONSTANTS(:,52).*(CONSTANTS(:,58)+STATES(:,28)./CONSTANTS(:,48))+STATES(:,29)))+( CONSTANTS(:,53).*CONSTANTS(:,57).*STATES(:,20))./( CONSTANTS(:,54).*(CONSTANTS(:,58)+STATES(:,23)./CONSTANTS(:,55)+STATES(:,24)./CONSTANTS(:,56)+STATES(:,29)./CONSTANTS(:,50))+STATES(:,20));
    RATES(:,20) = ( CONSTANTS(:,51).*STATES(:,29))./( CONSTANTS(:,52).*(CONSTANTS(:,58)+STATES(:,28)./CONSTANTS(:,48))+STATES(:,29)) - ( CONSTANTS(:,53).*CONSTANTS(:,57).*STATES(:,20))./( CONSTANTS(:,54).*(CONSTANTS(:,58)+STATES(:,23)./CONSTANTS(:,55)+STATES(:,24)./CONSTANTS(:,56)+STATES(:,29)./CONSTANTS(:,50))+STATES(:,20));
    RATES(:,30) =  - (( CONSTANTS(:,60).*STATES(:,21).*STATES(:,30))./( CONSTANTS(:,61).*(CONSTANTS(:,70)+STATES(:,23)./CONSTANTS(:,65))+STATES(:,30)))+( CONSTANTS(:,62).*CONSTANTS(:,59).*STATES(:,23))./( CONSTANTS(:,63).*(CONSTANTS(:,70)+STATES(:,24)./CONSTANTS(:,67)+STATES(:,29)./CONSTANTS(:,68)+STATES(:,20)./CONSTANTS(:,69))+STATES(:,23));
    RATES(:,23) = ((( CONSTANTS(:,60).*STATES(:,21).*STATES(:,30))./( CONSTANTS(:,61).*(CONSTANTS(:,70)+STATES(:,23)./CONSTANTS(:,65))+STATES(:,30)) - ( CONSTANTS(:,62).*CONSTANTS(:,59).*STATES(:,23))./( CONSTANTS(:,63).*(CONSTANTS(:,70)+STATES(:,24)./CONSTANTS(:,67)+STATES(:,29)./CONSTANTS(:,68)+STATES(:,20)./CONSTANTS(:,69))+STATES(:,23))) - ( CONSTANTS(:,64).*STATES(:,21).*STATES(:,23))./( CONSTANTS(:,65).*(CONSTANTS(:,70)+STATES(:,30)./CONSTANTS(:,61))+STATES(:,23)))+( CONSTANTS(:,66).*CONSTANTS(:,59).*STATES(:,24))./( CONSTANTS(:,67).*(CONSTANTS(:,70)+STATES(:,23)./CONSTANTS(:,63)+STATES(:,29)./CONSTANTS(:,68)+STATES(:,20)./CONSTANTS(:,69))+STATES(:,24));
    RATES(:,24) = ( CONSTANTS(:,64).*STATES(:,21).*STATES(:,23))./( CONSTANTS(:,65).*(CONSTANTS(:,70)+STATES(:,30)./CONSTANTS(:,61))+STATES(:,23)) - ( CONSTANTS(:,66).*CONSTANTS(:,59).*STATES(:,24))./( CONSTANTS(:,67).*(CONSTANTS(:,70)+STATES(:,23)./CONSTANTS(:,63)+STATES(:,29)./CONSTANTS(:,68)+STATES(:,20)./CONSTANTS(:,69))+STATES(:,24));
    RATES(:,31) =  - (( CONSTANTS(:,72).*STATES(:,24).*STATES(:,31))./( CONSTANTS(:,73).*(CONSTANTS(:,80)+STATES(:,32)./CONSTANTS(:,77))+STATES(:,31)))+( CONSTANTS(:,74).*CONSTANTS(:,71).*STATES(:,32))./( CONSTANTS(:,75).*(CONSTANTS(:,80)+STATES(:,33)./CONSTANTS(:,79))+STATES(:,32));
    RATES(:,33) = ( CONSTANTS(:,76).*STATES(:,24).*STATES(:,32))./( CONSTANTS(:,77).*(CONSTANTS(:,80)+STATES(:,31)./CONSTANTS(:,73))+STATES(:,32)) - ( CONSTANTS(:,78).*CONSTANTS(:,71).*STATES(:,33))./( CONSTANTS(:,79).*(CONSTANTS(:,80)+STATES(:,32)./CONSTANTS(:,75))+STATES(:,33));
    RATES(:,32) = ((( CONSTANTS(:,72).*STATES(:,24).*STATES(:,31))./( CONSTANTS(:,73).*(CONSTANTS(:,80)+STATES(:,32)./CONSTANTS(:,77))+STATES(:,31)) - ( CONSTANTS(:,74).*CONSTANTS(:,71).*STATES(:,32))./( CONSTANTS(:,75).*(CONSTANTS(:,80)+STATES(:,33)./CONSTANTS(:,79))+STATES(:,32))) - ( CONSTANTS(:,76).*STATES(:,24).*STATES(:,32))./( CONSTANTS(:,77).*(CONSTANTS(:,80)+STATES(:,31)./CONSTANTS(:,73))+STATES(:,32)))+( CONSTANTS(:,78).*CONSTANTS(:,71).*STATES(:,33))./( CONSTANTS(:,79).*(CONSTANTS(:,80)+STATES(:,32)./CONSTANTS(:,75))+STATES(:,33));
   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
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