N = 100;

dieAfter = 10;
initialInfection = 0.001;
infectionRate=0.4;
healRate=0.1;
birthRate=0.1;
%everyone on the boundary is immune
pop = zeros(N); 

%for starters, set everyone to healthy, and infect a few
pop(1:N,1:N) = ones(N) + (rand(N) < initialInfection); 

%create a temporary array
newpop  = pop;

%create colourmap
%you can be immune, healty or at different stages of the disease (up to
%being dead)
map = zeros(dieAfter+2,3);
map(1,:) = [0 1 0]; % zero is immune, colour green
map(2,:) = [1 1 1]; % one is healthy, colour white
map(3:end, 1) = flipud(linspace(1,0,size(map(3:end, 1),1))); % shades of red (up to black)

for t= 1:100
    t
    %show current state
    imagesc(pop,[0 dieAfter+2]);
    colormap(map);
   for i=1:N
       for j=1:N
          guy=pop(i, j);
          newguy=newpop(i, j);
          if guy>=dieAfter
              
              egeszs=0;
              for i2=-1:1 
                  
                       if i==1 && i2==-1
                          i2b=i2+N;
                      elseif i==N && i2==1
                          i2b=i2-N;
                       else i2b=i2;
                      end 

                  
                  for j2=-1:1
                      
                      if j==1 && j2==-1
                          j2b=j2+N;
                      elseif j==N && j2==1
                          j2b=j2-N;
                          else j2b=j2;
                      end 
                      
                      if pop(i+i2b, j+j2b)==0 || pop(i+i2b, j+j2b)==1
                          egeszs=egeszs+1;
                         if rand>birthRate && egeszs>=2;
                             newguy=1;
                         end
                      end                      
                  end
              end
              
          elseif guy==0
          elseif guy==1 %fertheto
              korny=0;
              for i1=-1:1 
                  
                       if i==1 && i1==-1
                          i1b=i1+N;
                      elseif i==N && i1==1
                          i1b=i1-N;
                          else i1b=i1;
                      end 
                  
                  for j1=-1:1
                      
                      if j==1 && j1==-1
                          j1b=j1+N;
                      elseif j==N && j1==1
                          j1b=j1-N;
                          else j1b=j1;
                      end 
                      
                      if pop(i+i1b, j+j1b)>1 && pop(i+i1b, j+j1b)<dieAfter
                          korny=korny+1;
                      end
                  end
              end
              if korny*rand>infectionRate
                  newguy=2;
              end
              
              %infectionRate
              
          elseif guy>1 %beteg
              if rand<healRate%%beall
                  newguy=0;
              else
                  newguy=guy+1;
              end
          end
          newpop(i,j)=newguy;
       end
   end
    
    
    pause(1);
    %update population with the temporary array
    pop = newpop;
    
end