Graph generated from Mathematica

 Clear[IC, a, b, h, i, m, n, f, NN, β, γ];
 IC = {499, 1}; 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_] := -((β I S)/NN) + γ I;
 f[2, t_, S_, I_] := (β I S)/NN - γ 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];
 Maroon = RGBColor[.5, 0, 0];
 ListPlot[t1, PlotMarkers -> {Automatic, Scaled[0.01]}, PlotStyle -> {DkYellow, Maroon}, AspectRatio -> 420/560]