function [S,I,R,Tout] = SIR_deterministic(beta,gamma,N,i_0,t_max) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the SIR model differential equations 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Inputs are:
% beta - infection rate, scalar
% gamma - recovery rate, scalar
% N - population size, scalar
% i_0 - initial number of infectives, scalar
% t_max - temporal horizon, scalar

% Outputs are:
% Tout - time instants in the interval [0,t_max], over which the solutions
%        are computed, vector with length specified automatically  
% S - evolution of susceptibles in [0,t_max], length(Tout)x1 vector
% I - evolution of infectives in [0,t_max], length(Tout)x1 vector
% R - evolution of removed in [0,t_max], length(Tout)x1 vector

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The state is Y = [S,I]'

% Define the differential equation right-hand side
SIRdiffeq = @(t, Y) [-beta*Y(1)*Y(2);
                    beta*Y(1)*Y(2) - gamma*Y(2)];

%% Initial condition
Y0 = [N-i_0,i_0]';

% Solve the differential equation (if needed check help ode45)
[Tout, Yout] = ode45(SIRdiffeq, [0 t_max], Y0);

% Return the desired functions S(t), I(t), R(t)
S = Yout(:,1);
I = Yout(:,2);
R = N - S - I;
end