There are a few builtin fcn in matlab to solve such ODEs.
% Symbolic
syms x(t) y(t) z(t) p k a1 a2 b2 c d sigma beta
Eqn1 = diff(x(t), t) == p*x*(1-y/k)-a1*x*y;
Eqn2 = diff(y(t), t) == c*a1*x*y-d*y-a2*y*z/(y+a2);
Eqn3 = diff(z(t), t) == sigma*z^2-beta*z^2/(y+b2);
Sols = dsolve([Eqn1, Eqn2, Eqn3])
%% Numerical solutions
ICs = [pi; pi/2; 2]; % Initial Conditions
[time, Sol]=ode45(@(t, xyz)FCN(t, xyz), [0, 10], ICs);
plot(time, Sol(:,1), 'r',time, Sol(:,2), 'g', time, Sol(:,3), 'b', 'linewidth', 2)
legend('x(t)', 'y(t)', 'z(t)'); grid on
title('Numerical Solutions')
xlabel('$t$', 'interpreter', 'latex')
ylabel('$x(t), \ y(t), \ z(t)$', 'interpreter', 'latex')
function dxyz = FCN(t, xyz)
p =.1;
k =.2;
a1 =.3;
a2 =.4;
b2 =.5;
c =.6;
d =.7;
sigma =.8;
beta =.9;
dxyz=[p*xyz(1)*(1-xyz(2)/k)-a1*xyz(1).*xyz(2);
c*a1*xyz(1).*xyz(2)-d*xyz(2)-a2*xyz(2).*xyz(3)./(xyz(2)+a2);
sigma*xyz(3).^2-(beta*xyz(3).^2)./(xyz(2)+b2)];
end
