Why not just solve them numerically:
% Parameters
tau = 0.138; % dimensionless kinetic parameter
delta = 9.9; % dimensionless kinetic parameter
phi = -0.456; % dimensionless kinetic parameter
% Initial conditions
X(1) = 0.178333333; % initial biomass concentration in g/L from experimental data
N(1) = 35.01666667; % initial culture nitrate concentration in mg/L from experimental data
q(1) = 79.85; % initial nitrogen quota in mg/g
f(1) = 120.3100949; % initial FAME yield in %wt from experimental data
B0 = [X(1) N(1) q(1) f(1)];
tspan = [0 200];
[t, B] = ode45(@odefn, tspan, B0);
X = B(:,1); N = B(:,2); q = B(:,3); f = B(:,4);
% Chlorophyll fluorescence
Y_II = exp(tau*q)./(exp(tau*q)+delta) + phi;
subplot(2,1,1)
plot(t,X),grid
xlabel('t'),ylabel('Algal biomass growth rate')
subplot(2,1,2)
plot(t,N),grid
xlabel('t'),ylabel('Nitrate consumption')
figure
subplot(2,1,1)
plot(t,q),grid
xlabel('t'),ylabel('Accum. rate of nitrogen quota')
subplot(2,1,2)
plot(t,f),grid
xlabel('t'),ylabel('FAME production')
figure
plot(t,Y_II),grid
xlabel('t'),ylabel('Chlorophyll fluorescence')
function dBdt = odefn(~,B)
% Data
mu_d = 0; % specific decay rate
mu_N = 2.69; % maximum specific nitrate uptake rate
K_N = 0.8; % half velocity coefficient
k_q = 19.6; % minimum nitrogen quota
theta = 6.69; % modified kinetic constant for biolipid syntheiss
epsilon = 0.001; % modified parameter taking into account the complex effect of lipid esterification
gamma = 7.53; % modified kinetic constant for biolipid consumption
I_0 = 80; % incident light intensity micromol/m^2.s from experimental data
alpha = 196.4; % cell absorption coefficient
beta = 0; % bubble scattering coefficient
L = 0.044; % reactor length
k_i = 100; % light inhibition terms
mu_max = 0.36; % maximum specific growth rate in 1 h
k_s = 91.2; % light saturation term in micromol/m^2.s
% Extract parameters
X = B(1); N = B(2); q = B(3); f = B(4);
% Modified Beer-Lambert law
I = zeros(11,1);
for z = 0:10
I(z+1) = I_0 * exp(-(alpha*X + beta) * (z*L)/10);
end
% Trapezoidal rule
mu_m_first = I(1)/(I(1)+k_s+I(1)^2/k_i);
mu_m_last = I(11)/(I(11)+k_s+I(11)^2/k_i);
mu_m_middle = 0;
for n = 2:10
mu_m_middle = mu_m_middle + (I(2)/(I(2)+k_s+I(2)^2/k_i));
end
mu_m = mu_max/20 * (mu_m_first + 2*mu_m_middle + mu_m_last);
% Specific growth rate
mu_0 = mu_m * (1 - k_q/q);
% Rate equations
dBdt = [mu_0*X - mu_d*X;
-mu_N * (N/(N+K_N)) * X;
mu_N * (N/(N+K_N)) - mu_m * (1-k_q/q) * q;
mu_m * (theta*q - epsilon*f) * (1-k_q/q) - gamma*mu_N*(N/(N+K_N))];
end
