function[A] = toom(L,T,p)
A = zeros(T,2*L);
% Permute the numbers {1,...,L}
% This is fractionally easier than changin the colour map
list = randperm(L);
for i = 1:L
    % Each site is initialised with a different type
    A(T,2*i) = list(i);
end

for t=T-1:-1:1
    for i=1:L
        % Position in matrix to update
        x = 2*i - mod(t,2);
        % Left parent of x
        left = mod(x-2,2*L) + 1; 
        % Right parent of x
        right = mod(x,2*L) + 1;
        % Check parents type
        if A(t+1,left) == A(t+1,right)
            A(t,x) = A(t+1,left);
        else
            if rand > p
                A(t,x) = A(t+1,left);
            else
                A(t,x) = A(t+1,right);
            end
        end
    end
end

% The following is only to make the image look nice
colordata = colormap;
colordata(1,:) = [0 0 0]; % Lowest number gets mapped to black
colormap(colordata);
%--------------------------------------------------

imagesc(A)