function GOL(n,m,tmax)

A = sparse(round(rand(n,m)));
B = sparse(size(A));
spy(A); pause(0.05)

for t = 1:tmax
    
    for i = 1:n %row index
        for j = 1:m %column index
            %find the cell's neighbours
            try
                %same row, to the left and right of cell
                %row above cell, to the left, direcly above and to the right
                %row below cell, to the left, direcly below and to the right
                neighbours = [A(i,j-1), A(i,j+1), A(i-1, j-1), A(i-1, j), A(i-1, j+1),A(i+1, j-1), A(i+1, j), A(i+1, j+1)];
            catch
                %indecies are out of bounds
                %apply periodic boundary conditions:
                % the top and bottom row are neighbours (fold matrix into a tube)
                % the left and right edges of the matrix are neighbours (make the tube into a doughnut)
                if i == 1
                    if j ==1
                        neighbours = [A(i,m), A(i,j+1), A(n, m), A(n, j), A(n, j+1),A(i+1, m), A(i+1, j), A(i+1, j+1)];
                    elseif j==m
                        neighbours = [A(i,j-1), A(i,1), A(n, j-1), A(n, j), A(n, 1),A(i+1, j-1), A(i+1, j), A(i+1, 1)];
                    else
                        neighbours = [A(i,j-1), A(i,j+1), A(n, j-1), A(n, j), A(n, j+1),A(i+1, j-1), A(i+1, j), A(i+1, j+1)];
                    end
                elseif i == n
                    if j ==1
                        neighbours = [A(i,m), A(i,j+1), A(i-1, m), A(i-1, j), A(i-1, j+1),A(1, m), A(1, j), A(1, j+1)];
                    elseif j == m
                        neighbours = [A(i,j-1), A(i,1), A(i-1, j-1), A(i-1, j), A(i-1, 1),A(1, j-1), A(1, j), A(1, 1)];
                    else
                        neighbours = [A(i,j-1), A(i,j+1), A(i-1, j-1), A(i-1, j), A(i-1, j+1),A(1, j-1), A(1, j), A(1, j+1)];
                    end
                else
                    if j ==1
                        neighbours = [A(i,m), A(i,j+1), A(i-1, m), A(i-1, j), A(i-1, j+1),A(i+1, m), A(i+1, j), A(i+1, j+1)];
                    elseif j==m
                        neighbours = [A(i,j-1), A(i,1), A(i-1, j-1), A(i-1, j), A(i-1, 1),A(i+1, j-1), A(i+1, j), A(i+1, 1)];
                    else
                        display (['I dont know what the error is 3 :( i=' num2str(i) ' j=' num2str(j)]) % the error was not due to out of bounds indecies
                        return; %quit gracefully
                    end
                end
            end
            
            %count live neighbours
            live = 0;
            for k = 1:length(neighbours)
                if neighbours(k)==1
                    live = live+1;
                end
            end
            
            if A(i,j)==0 %the cell is dead
                %if exactly 3 live neighbours, the cell becomes alive
                if live == 3
                    B(i,j) = 1;
                end
            else % the cell is alive
                %if less then 2 live neighbours, cell dies
                if live < 2
                    B(i,j) = 0;
                %if more then 3 live neighbours, cell dies
                elseif live > 3
                    B(i,j) = 0;
                %otherwise it lives
                else
                    B(i,j) = 1;
                end                  
            end
        end
    end
    
    A = B;
    spy(A); pause(0.05)
end