Summary

Clear[IC,a,b,h,i,m,n,f,NN,β,γ];
IC={3000000,1,0};m=Length[IC];
a=0;b=60;h=0.5;
NN=1;
n=(b-a)/h;
β=0.001;γ=0.1;
f[1,t_,S_,I_,R_]:=-((βIS)/NN);
f[2,t_,S_,I_,R_]:=(βIS)/NN-γI;
f[3,t_,S_,I_,R_]:=γI;
Clear[t,k1,k2,k3,k4,w];
t[i_]:=a+i*h;
For[i=1,i<=m,i++,w[i,0]=IC[[i]]];
k1[i_,j_]:=h*(f@@Flatten[{i,t[j-1],Table[w[k,j-1],{k,1,m}]}]);
k2[i_,j_]:=h*(f@@Flatten[{i,t[j-1]+h/2,Table[w[k,j-1]+k1[k,j]/2,{k,1,m}]}]);
k3[i_,j_]:=h*(f@@Flatten[{i,t[j-1]+h/2,Table[w[k,j-1]+k2[k,j]/2,{k,1,m}]}]);
k4[i_,j_]:=h*(f@@Flatten[{i,t[j-1]+h,Table[w[k,j-1]+k3[k,j],{k,1,m}]}]);
w[i_,j_]:=w[i,j]=w[i,j-1]+(k1[i,j]+2*k2[i,j]+2*k3[i,j]+k4[i,j])/6
t1=Table[Table[{t[j],w[i,j]},{j,0,n}],{i,1,m}](*//TableForm*);
DkYellow=RGBColor[.808,.808,.192];
Teal=RGBColor[0,.5,.5];
Maroon=RGBColor[.5,0,0];
ListPlot[t1,PlotMarkers->{Automatic,Scaled[0.01]},PlotStyle->{DkYellow,Maroon,Teal},AspectRatio->420/560]

Licensing