function [S,I,R,t]=sir_solve(beta,gamma,t_offset,I_init)

% Parameters used to generate the simulation:
% S_init=4.9e6
% I_init=10
% beta=1
% gamma=1/5
% T_fin=100;
% t_offset=0;

S_init=4.9e6;
R_init=0;
T_fin=100;

N=S_init+I_init+R_init;

% options = odeset('RelTol',1e-8,'AbsTol',1e-8,'Stats','on','OutputFcn',@odeplot);
options = odeset('RelTol',1e-8,'AbsTol',1e-8);
[T,Y]=ode15s(@myfun,[t_offset,T_fin],[S_init;I_init;R_init],options);

S=Y(:,1);
I=Y(:,2);
R=Y(:,3);
t=T;

ti=0:1:T_fin;
Si=interp1(t,S,ti,'spline');
Ii=interp1(t,I,ti,'spline');
Ri=interp1(t,R,ti,'spline');

t=ti;
S=Si;
I=Ii;
R=Ri;
    
    function dydtau=myfun(~,y)
        dsdtau=-(beta/N)*y(1)*y(2);
        didtau= (beta/N)*y(1)*y(2)-gamma*y(2);
        drdtau= gamma*y(2);
    
        dydtau=[dsdtau;didtau;drdtau];
end

end