function f= fun2(~,x)
sigma1 = 0.6; sigma2 = 5; sigma3 = 0.5; sigma4 = 2;
r1 = 1; r2 = 0.5; r3 = 0.6;
N1 = 1000; N2 = 300; N3 = 20;
f = @f;
f = [r1*x(1)*(1-x(1)/N1-sigma1*x(2)/N2);
r2*x(2)*(-1-x(2)/N2 + sigma2*x(1)/N1-sigma3*x(3)/N3);
r3*x(3)*(-1-x(3)/N3 + sigma4*x(2)/N2)];
[t,x] = ode45 ('fun2', [0, 15], [900, 200, 15]);
plot(t,x), grid