Customizable S-i-R code

This code lets you set S, i, R, starting time, ending time, and number of increments within the time interval.  It prints out t, S, i, R for every time-step.

Clear[t];
Clear[tinitial];
Clear[tfinal];
Clear[numberofsteps];
Clear[S];
Clear[i];
Clear[R];
Clear[deltat];
S =45400;
i = 2100;
R = 2500;
tinitial=0;
tfinal = 1;
t = tinitial;
numberofsteps = 10;
deltat = (tfinal- tinitial)/numberofsteps;
For[k=1,k≤numberofsteps,
Sprime = -.00001*S*i; 
iprime =.00001*S*i - i/14;
 Rprime = i/14; 
deltaS = Sprime*deltat; 
deltai = iprime*deltat; 
deltaR = Rprime*deltat; 
t=t+deltat; 
S = S+deltaS; 
i = i + deltai; 
R = R+deltaR; 
Print[N[t],"   ", S,"    ",i,"    ", R];
 k++]

This code is identical to the above, except it only prints out values of S, i, R for the ending time.

Clear[t];
Clear[tinitial];
Clear[tfinal];
Clear[numberofsteps];
Clear[S];
Clear[i];
Clear[R];
Clear[deltat];
S =45400;
i = 2100;
R = 2500;
tinitial=0;
tfinal = 1;
t = tinitial;
numberofsteps = 10;
deltat = (tfinal- tinitial)/numberofsteps;
For[k=1,k≤numberofsteps,
Sprime = -.00001*S*i; 
iprime =.00001*S*i - i/14;
 Rprime = i/14; 
deltaS = Sprime*deltat; 
deltai = iprime*deltat;
 deltaR = Rprime*deltat; 
t=t+deltat; 
S = S+deltaS; 
i = i + deltai; 
R = R+deltaR; 
 k++]
Print[t];Print[S];Print[i];Print[ R];

There and Back

Suppose you want to measure the accumulated error from doing Euler's method (remember, that's the official name for our S-i-R calculations).  What you want to have the computer do is 
(1) run the S-i-R code from tinitial to tfinal, 
(2) note the ending S-i-R from that procedure 
(3) use them as the initial values for running the S-i-R code from old tfinal to old tinitial.  Here's one way to do it:

Clear[t];
Clear[tinitial];
Clear[tfinal];
Clear[numberofsteps];
Clear[S];
Clear[i];
Clear[R];
Clear[deltat];
S =45400;
i = 2100;
R = 2500;
tinitial=0;
tfinal = 1;
t = tinitial;
numberofsteps = 10;
deltat = (tfinal- tinitial)/numberofsteps;
For[k=1,k≤numberofsteps,
Sprime = -.00001*S*i; 
iprime =.00001*S*i - i/14;
 Rprime = i/14; 
deltaS = Sprime*deltat; 
deltai = iprime*deltat;
 deltaR = Rprime*deltat; 
t=t+deltat; 
S = S+deltaS; 
i = i + deltai; 
R = R+deltaR; 
 k++]
Print[t];Print[S];Print[i];Print[ R];

Clear[t];
Clear[tinitial];
Clear[tfinal];
Clear[numberofsteps];
Clear[S];
Clear[i];
Clear[R];
Clear[deltat];
S =put the S output here;
i = put the i output here;
R = put the R output here;
tinitial=1;
tfinal = 0;
t = tinitial;
numberofsteps = 10;
deltat = (tfinal- tinitial)/numberofsteps;
For[k=1,k≤numberofsteps,
Sprime = -.00001*S*i; 
iprime =.00001*S*i - i/14;
 Rprime = i/14; 
deltaS = Sprime*deltat; 
deltai = iprime*deltat;
 deltaR = Rprime*deltat; 
t=t+deltat; 
S = S+deltaS; 
i = i + deltai; 
R = R+deltaR; 
 k++]
Print[t];Print[S];Print[i];Print[ R];

Oh, but there's a more efficient way to do it.  The thing is that you can run it without re-setting S, i, and R!  See below.

Clear[t];
Clear[tinitial];
Clear[tfinal];
Clear[numberofsteps];
Clear[S];
Clear[i];
Clear[R];
Clear[deltat];
S =45400;
i = 2100;
R = 2500;
tinitial=0;
tfinal = 1;
t = tinitial;
numberofsteps = 10;
deltat = (tfinal- tinitial)/numberofsteps;
For[k=1,k≤numberofsteps,
Sprime = -.00001*S*i; 
iprime =.00001*S*i - i/14;
 Rprime = i/14; 
deltaS = Sprime*deltat; 
deltai = iprime*deltat;
 deltaR = Rprime*deltat; 
t=t+deltat; 
S = S+deltaS; 
i = i + deltai; 
R = R+deltaR; 
 k++]
tinitial=1;
tfinal = 0;
t = tinitial;
deltat = (tfinal- tinitial)/numberofsteps;
For[k=1,k≤numberofsteps,
Sprime = -.00001*S*i; 
iprime =.00001*S*i - i/14;
 Rprime = i/14; 
deltaS = Sprime*deltat; 
deltai = iprime*deltat;
 deltaR = Rprime*deltat; 
t=t+deltat; 
S = S+deltaS; 
i = i + deltai; 
R = R+deltaR; 
 k++]
Print[t];Print[S];Print[i];Print[ R];

Plotting S-i-R data (just the basics)

This version of the code is identical to the code in the previous section, except that it takes all of the values for t, S, i, R and puts them in a list so that we can plot them later.

Clear[t];
Clear[tinitial];
Clear[tfinal];
Clear[numberofsteps];
Clear[S];
Clear[i];
Clear[R];
Clear[deltat];
Clear[datapoints];
datapoints = {};
S =45400;
i = 2100;
R = 2500;
tinitial=0;
tfinal = 1;
t = tinitial;
numberofsteps = 10;
deltat = (tfinal- tinitial)/numberofsteps;
For[k=1,k≤numberofsteps,
Sprime = -.00001*S*i; 
iprime =.00001*S*i - i/14; 
Rprime = i/14; 
deltaS = Sprime*deltat;
 deltai = iprime*deltat;
 deltaR = Rprime*deltat;
 t=t+deltat; 
S = S+deltaS; 
i = i + deltai; 
R = R+deltaR;
datapoints =Append[datapoints, {t,S,i,R}]; k++]

Notice that there is no output when you run this code.  That's so that you don't get a huge long mess in cases where numberofsteps is large.  But if you want to see the list you generated, just say...

datapoints

And of course you can change any of the variable names, as long as you change each occurance of the name.  You might want to do this if you have more than one set of data you'd like to generate, so that you don't overwrite your old data when you run the program again.

Here's how we make a plot of S from the data:

tversusS = Table[{datapoints[[j,1]],datapoints[[j,2]]}, {j,1, Length[datapoints]}];
ListPlot[tversusS]

The double-brackets tell Mathematica that you're looking inside the list; the j is a counter which goes through each chunk of the list (each part that's set off in {}), and the 1 and 2 refer to the first and second slots in each {t, S, i, R} chunk.  To plot the infecteds, substitute 3 for 2, and to plot the recovereds, substitute 4 for 2 in the above command.  
Oh, and here's how to (a) do it all in one step and (b) connect the dots:

ListPlot[Table[{datapoints[[i,1]],datapoints[[i,2]]},{i,1,Length[datapoints]}],Joined->True]

Did you want to label your plot?  Here's how:

ListPlot[tversusS, PlotLabel -> "Crazy!"]

...well, except that you probably want a different label.
Tune in to the next notebook to find out more interesting things Mathematica will let you do with plotting.