N=1000; I0=10; beta=0.1; gamma=0.05; tend=300;
tspan = [0,tend];

S0 = N - I0;

R0 = 0;

y0 = [S0; I0; R0];

[t,y] = ode45(@(t,y) SIRfunc(t,y,beta,gamma,N), tspan, y0);

plot(t,y(:,1),'-',t,y(:,2),'m-',t,y(:,3),'g-', 'LineWidth' , 2);

legend('Susceptible','Infectious','Recovered')

title('SIR model')

function dydt = SIRfunc(t,y,beta,gamma,N)

dydt = [-beta/N*y(2)*y(1); ...

beta/N*y(2)*y(1) - gamma*y(2); ...

gamma*y(2)];

end