nhealthy = 5*10^3; nill = 10; nimmune = 0; max = nhealthy;
illradius = 0.1; step = 0.1;
healthy0 = RandomReal[{-10, 10}, {nhealthy, 2}]; ill0 = RandomReal[{-10, 10}, {nill, 2}]; immune0 = {};
healthy = healthy0; ill = ill0; immune = immune0;
probilldie = 0.01; probinfect = 0.25; probheal = 0.005;
evo = Reap[
    (*Who gets ill*)
    Do[
     For[j = 1, j <= nill, j++,
      infections = Position[Norm[# - ill[[j]]] & /@ healthy, x_ /; (x < illradius)]; (*Potential infections*)
      random = RandomReal[{0, 1}, Dimensions[infections][[1]] ];
      infections = Position[random, x_ /; (x < probinfect)];(*actual infections*)
      ill = Join[ill, Extract[healthy, infections]];
      healthy = Delete[healthy, infections]; (*moving the infected from the healthy to the ill*)
      ];
     nhealthy = Dimensions[healthy][[1]]; 
     nill = Dimensions[ill][[1]]; nimmune = Dimensions[immune][[1]];
     
     (*Who heals and get immune*)
     random = RandomReal[{0, 1}, Dimensions[ill][[1]] ];
     healings = Position[random, x_ /; (x < probheal)];
     immune = Join[immune, Extract[ill, healings] ];
     ill = Delete[ill, healings];
     nhealthy = Dimensions[healthy][[1]]; 
     nill = Dimensions[ill][[1]]; nimmune = Dimensions[immune][[1]];
     
     (*Who dies*)
     random = RandomReal[{0, 1}, Dimensions[ill][[1]] ];
     deaths = Position[random, x_ /; (x < probilldie)];
     ill = Delete[ill, deaths];
     nhealthy = Dimensions[healthy][[1]]; 
     nill = Dimensions[ill][[1]]; nimmune = Dimensions[immune][[1]];
     
     (*Movement*)
     healthy = Clip[healthy + RandomVariate[NormalDistribution[0, step], {nhealthy, 2}], {-10, 10}];
     ill = Clip[ill + RandomVariate[NormalDistribution[0, step], {nill, 2}], {-10, 10}];
     immune = Clip[immune + RandomVariate[NormalDistribution[0, step], {nimmune, 2}], {-10, 10}];
     
     Sow[{healthy, ill, immune}];
     If[nill < 1, Break[]];
     , {500}]
    ][[2, 1]];
frames = Table[
   Grid[{{
      Graphics[{PointSize[0.01], Point[evo[[j, 1]]], Red, Thick, Circle[#, illradius] & /@ evo[[j, 2]] , Green, Point[evo[[j, 3]]]}, PlotRange -> {{-10, 10}, {-10, 10}}, ImageSize -> Medium]
      ,
      ListPlot[{(Dimensions /@ evo[[All, 1]])[[1 ;; j, 1]], (Dimensions /@ evo[[All, 2]])[[1 ;; j, 1]], (Dimensions /@ evo[[All, 3]])[[1 ;; j, 1]], (Dimensions /@ evo[[All, 1]])[[1 ;; j, 1]] + (Dimensions /@ evo[[All, 2]])[[1 ;; j, 1]] + (Dimensions /@ evo[[All, 3]])[[1 ;; j, 1]]}, PlotRange -> {{0, Dimensions[evo][[1]]}, {0, max*1.1}}, PlotStyle -> {Black, Red, Green, Blue}, Axes -> False, Frame -> True, FrameTicks -> None, FrameLabel -> {"Time", "#"}, LabelStyle -> {Black, Bold}, Joined -> True, PlotLegends -> LineLegend[{Black, Red, Green, Blue}, {"Susceptible", "Ill", "Immune", "Total population"}], ImageSize -> Medium ]
      }}]
   , {j, 1, Dimensions[evo][[1]], 2}];
ListAnimate[frames]