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_] := -((β I S)/NN);
   f[2, t_, S_, I_, R_] := (β I S)/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]