nuy = nuy_max*S/(Ks + S);
eqn2 = diff(X) == (nuy - Q/V)*X;
eqn3 = diff(P*V) == nuy_g*X*V - K*P*V;
eqn4 = diff(S*V) == S_F*Q - nuy*X*V/Y;
eqns = [eqn1 eqn2 eqn3 eqn4];
[V, S] = odeToVectorField(eqns)
V =

S =

odefun = matlabFunction(V, 'vars', {'t','Y'});
[t, y] = ode45(odefun, tspan, y0);
tb = tiledlayout(2, 2, 'TileSpacing', 'Compact');
plot(t, y(:,1)), grid on, title('X(t)')
plot(t, y(:,2)), grid on, title('V(t)')
plot(t, y(:,3)), grid on, title('P(t)')
plot(t, y(:,4)), grid on, title('S(t)')