ode.m

% General constants obtained from Section 3.1
mu_X_max    = 0, 275
K_XS        = 0.0624
K_iXS       = 12.40
K_XN        = 0.0812
K_iXN       = 0.6250
K_I         = 45.50
K_iI        = 297.75
A_0X        = 0.6175
E_aX        = 25.9243
B_0X        = 0.1101
E_bX        = 48.0151
q_L_max      = 0.1512
K_LS        = 4.9155
K_iLS       = 0.1375
K_iNL       = 380.31
A_0L        = 0.7
E_aL        = 20
B_0L        = 0.08
E_bL        = 32
Y_XS        = 1.1025
Y_LS        = 0.048
Y_XN        = 5.1623
K_h         = 0.7606
rho         = 0.7419
l           = 0.25
T_0         = 293

function output = f_X(T)
    % Equation 6
    output = A_0X * exp( (-E_aX / R) * (1/T - 1/T_0) ) - B_0X * exp( (-E_bX / R) * (1/T - 1/T_0) )
end

function output = f_L(T)
    % Equation 7
    output = A_0L * exp( (-E_aL / R) * (1/T - 1/T_0) ) - B_0L * exp( (-E_bL / R) * (1/T - 1/T_0) )
end

function diffeqs = ode_sys(t, var)
    % t [vector]: range or array of time values 
    % var [vector]: array of input variables

    %% Extract dependent variables from var
    X = var(1);     % Used in Equation 8, and 9
    L = var(2);     % Used in Equation 9
    S = var(3);     % Used in Equation 10
    N = var(4);     % Used in Equation 11
    H = var(5);     % Used in Equation 12

    % Equation 4
    mu_X = mu_X_max * f_X(S) * f_X(N) * f_X(I) * f_X(T)

    % Equation 5
    mu_L = q_L_max * f_L(S) * f_L(N) * f_L(I) * f_L(T)

    % dX/dt [Equation 8]
    diffeqs(1, 1) = (mu_X * X);

    % dL/dt [Equation 9]
    diffeqs(1, 2) = (mu_L * X);

    % dS/dt [Equation 10]
    diffeqs(1, 3) = ((-1 / Y_XS) * diff(X, t)) - ((1 / Y_LS) * diff(L, t));
    
    % dN/dt [Equation 11]
    diffeqs(1, 4) = (-1/Y_XN) * diff(X, t)

    % dH/dt [Equation 12]
    diffeqs(1, 5) = -K_h * diff(S, t)
end

range = [1:50]
initial_conditions = [600, 3, 0, 3, 3]

[time_solution, var_solution] = ode45(@ode_sys, range, initial_conditions)

%figure
plot (time_solution, var_solution)
title('kinetic model for microalgae')
xlabel('time,day')
ylabel('growth,g')