%this function uses Euler's method to forwars/backward solve
%the SIR Model
function [S,I,R,t] = fwdbwd(tmax,N,a,b,S0,I0,R0)
%when calling this function, you need to pass the maximum time (tmax),
%the number of time steps (N), the parameters a and b, and the initial 
%populations S0,I0, and R0.
%note that delta t = tmax/N
delt = tmax/N; %this determines the size of the time step delta t
S = zeros(1,N+1); %define an empty susceptible vector
I = zeros(1,N+1); %define an empty infectives vector
R = zeros(1,N+1); %define an empty recovered vector
S(1) = S0; %first value in the Susceptible vector is the inital cond.
I(1) = I0; %first value in the Infectives vector is the inital cond.
R(1) = R0; %first value in the Recovered vector is the inital cond.
for n = 1:N %this is the iterative loop for Euler's Method
    S(n+1) = S(n) - a*S(n)*I(n)*delt;%D.E. for Susceptibles
    I(n+1) = I(n) + (a*S(n)*I(n) - b*I(n))*delt;%D.E. for Infectives
    R(n+1) = R(n) + b*I(n)*delt;%D.E. for Recovered
end %end the for loop
%Now we will step backwards
Sb = zeros(1,N+1); %define an empty susceptible vector for back stepping
Ib = zeros(1,N+1); %define an empty infectives vector for back stepping
Rb = zeros(1,N+1); %define an empty recovered vector for back stepping
Sb(1) = S(N+1); %Initialize Sb with S(tmax) value
Ib(1) = I(N+1); %Initialize Ib with I(tmax) value
Rb(1) = R(N+1); %Initialize Rb with R(tmax) value
delt = -delt;%Our delta t is now negative for back stepping
for n = 1:N %this is the iterative loop for Euler's Method
    Sb(n+1) = Sb(n) - a*Sb(n)*Ib(n)*delt;%D.E. for Susceptibles
    Ib(n+1) = Ib(n) + (a*Sb(n)*Ib(n) - b*Ib(n))*delt;%D.E. for Infectives
    Rb(n+1) = Rb(n) + b*Ib(n)*delt;%D.E. for Recovered
end %end the for loop

%we now create a table of our results
Initial = [S0;I0;R0];
BackValues = [Sb(N+1);Ib(N+1);Rb(N+1)];
rnames = {'S';'I';'R'};
T = table(Initial,BackValues,'RowNames',rnames)