Plotting the phase portrait for a system of 3 ODE by matplotlib

Pre-computing the solutions for the given equilibrium points as

sol=[ odeint(dynamical_model_bacteria,init,t) for init in x0];

you can then use them in separate loops to construct the two diagrams. Using subplot reduces the number of image windows, the assembled graphs then give

After that use the same axes object to draw all the curves in it, add some for some random initial points to fill the space

for k in range(len(x0)):
    x,y1,y2=sol[k].T
    ax.plot(x,y1,y2,'r', lw=2)

for k in range(80):
    sol = odeint(dynamical_model_bacteria,0.3*np.random.rand(3),t)
    x,y1,y2=sol.T
    ax.plot(x,y1,y2,'b', lw=1)

The saved plot then is

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.integrate import solve_ivp, odeint
s=0  
omega=np.array([2,3.7,4,3])
omega11,omega12,omega21,omega22=omega


beta=np.array([28,24])
beta1,beta2=beta

epsilon=np.array([1.12,3.1])
epsilon1,epsilon2=epsilon

g12=2
gamma12=2
def dynamical_model_bacteria(X,t):
    x, y1, y2 = X
    dxdt=x*(1-x-y1-y2)+s
    dy1dt=y1*(-epsilon[0]*(1+omega11*y1+omega12*y2)-g12*y2+beta1*x)
    dy2dt=y2*(-epsilon[1]*(1+omega21*y1+omega22*y2)+gamma12*y1+beta2*x) 
    return [dxdt,dy1dt,dy2dt]

print( dynamical_model_bacteria([50,50,50],0))  
# compute some solutions in "orderly" fashion
x0=np.array([
       [ 0.11111111,  0.88888889,  0.        ],
       [ 0.37237237,  0.        ,  0.62762763],
       [ 0.        ,  0.10813086, -0.22171438],
       [ 0.17247589,  0.35219856,  0.47532555]])
t=np.linspace(0,30,8000)

sol=[ odeint(dynamical_model_bacteria,init,t) for init in x0];
# first plot
fig=plt.figure(figsize=(15,5))
for k in range(len(x0)):
    print( dynamical_model_bacteria(x0[k],0)) 
    plt.subplot(2,2,k+1);
    x,y1,y2=sol[k].T
    plt.xlabel('Time  (day)')
    plt.ylabel('populations (num/ml)')
    plt.title(r' $\quad\beta_{1}<\epsilon_{1} and\ \beta_{2}<\epsilon_{2}$'+'      '+'$x_{0}=$' +str(x0[k]))
    plt.plot(t,x,'r--' ,linewidth=2.5, markersize=1.5)
    plt.plot(t,y1,'b:', linewidth=3, markersize=0.5)
    plt.plot(t,y2,'g-.', linewidth=3.5, markersize=0.2)
    # set the limits
    #plt.xlim([0, 30])
    plt.ylim([-1, 1])
    plt.legend( ('Bacteria','protozoa','Daphnia'),
           loc='upper center', shadow=True)
plt.tight_layout(); 
plt.savefig('/tmp/EP2_fig1.png', dpi=300, bbox_inches='tight')
#second plot, adding some random solutions with IC in the same size range
fig=plt.figure(figsize=(15,5))
ax = plt.axes(projection="3d")
for k in range(len(x0)):
    x,y1,y2=sol[k].T
    ax.plot(x,y1,y2,'r', lw=2)

for k in range(80):
    sol = odeint(dynamical_model_bacteria,0.3*np.random.rand(3),t)
    x,y1,y2=sol.T
    ax.plot(x,y1,y2,'b', lw=1)

plt.savefig('/tmp/EP2_fig2.png', dpi=300, bbox_inches='tight')
plt.show(); plt.close()