If you use my work, please let me know.  I would also love to see any changes or additions you might find necessary.
 

%A script to generate a patchy population.

%Initial inputs:
n  =  input('Number of iterations?');
g  =  input('Grid size?');
r  =  input('growth rate? (0-100%)')/100 + 1;
d  =  input('Rate of disturbance events? (0-100%)')/100;
c  =  input('Rate of colonization? (0-100%)')/100;
p  =  input('Patchiness factor? (0-100%)')/10;
t  =  input('topological complexity? (0-100%)')/100;

%Initial population on seafloor set to zero
seafloor = zeros(g,g);

%Seeding of organisms
for count = 1:(g/10)
    rndx = ceil(rand(1,1)*g);
    rndy = ceil(rand(1,1)*g);
    seafloor(rndx,rndy) = .1;
end

%To create geographical effects

    %Initial topography of seafloor set to zero
    g2 = ceil(g*t);
    topo = zeros(g2,g2);

    %initial hot spots
    for count = 1:(3*g2);
        rndx = ceil(rand(1,1)*g2);
        rndy = ceil(rand(1,1)*g2);
        topo(rndx,rndy) = .1;
    end

    %Begin iterative process
    for count = 1:n;

        %Deflation
        for x = 2:(g2-1);
            for y = 2:(g2-1);
                if topo(x,y) >= 10, topo(x,y)=0; topo(x+1,y)=0; topo(x+1,y+1)=0;
                                           topo(x,y+1)=0; topo(x-1,y)=0; topo(x-1,y-1)=0;
                                           topo(x,y-1)=0; topo(x+1,y-1)=0; topo(x-1,y+1)=0;
                end
            end
        end

        %Volcanism
        for x = 2:(g2-1);
            for y = 2:(g2-1);
                topo(x,y) = [topo(x+1,y)+topo(x+1,y+1)+topo(x,y+1)+
                                  topo(x-1,y)+topo(x-1,y-1)+topo(x,y-1)+
                                  topo(x+1,y-1)+topo(x-1,y+1)+topo(x,y)]/9*r;
            end
        end

end

%Convert 'topo' to a matrix of percentages
topo = topo/20;

%Begin population growth

    %Begin iterative process
    for count = 1:n

        %Die-off
        for x = 2:(g-1);
            for y = 2:(g-1);
                if seafloor(x,y) >= 10, seafloor(x,y)=seafloor(x,y)*topo(ceil(x*t),ceil(y*t));
                                                seafloor(x+1,y)=seafloor(x+1,y)*topo(ceil(t*(x+1)),ceil(t*y));
                                                seafloor(x+1,y+1)=seafloor(x+1,y+1)*topo(ceil(t*(x+1)),ceil(t*(y+1)));
                                                seafloor(x,y+1)=seafloor(x,y+1)*topo(ceil(t*x),ceil(t*(y+1)));
                                                seafloor(x-1,y)=seafloor(x-1,y)*topo(ceil(t*(x-1)),ceil(t*y));
                                                seafloor(x-1,y-1)=seafloor(x-1,y-1)*topo(ceil(t*(x-1)),ceil(t*(y-1)));
                                                seafloor(x,y-1)=seafloor(x,y-1)*topo(ceil(t*x),ceil(t*(y-1)));
                                                seafloor(x+1,y-1)=seafloor(x+1,y-1)*topo(ceil(t*(x+1)),ceil(t*(y-1)));
                                                seafloor(x-1,y+1)=seafloor(x-1,y+1)*topo(ceil(t*(x-1)),ceil(t*(y+1)));
                end
            end
        end

        %Disturbance algorithm
        for count = 1:(d*g)
            rndx = ceil(rand(1,1)*g);
            rndy = ceil(rand(1,1)*g);
            seafloor(rndx,rndy) = 0;
        end

        %Colonization algorithm
        for count = 1:(g*c)
            rndx = ceil(rand(1,1)*g);
            rndy = ceil(rand(1,1)*g);
            seafloor(rndx,rndy) = seafloor(rndx,rndy) + .1;
        end

        %Growth algorithm
        for x = 2:(g-1)
            for y = 2:(g-1)
                seafloor(x,y) = [seafloor(x+1,y)+seafloor(x+1,y+1)+seafloor(x,y+1)+seafloor(x-1,y)+
                                       seafloor(x-1,y-1)+seafloor(x,y-1)+seafloor(x+1,y-1)+seafloor(x-1,y+1)+
                                       seafloor(x,y)]/9*r;
            end
         end

    end

 %Sink the topography to create azoic zones
 for x = 1:g
     for y = 1:g
         if seafloor(x,y) <= p, seafloor(x,y) = 0;
         else seafloor(x,y) = seafloor(x,y) - p;
        end
    end
end

%Final plots
figure
contour(seafloor,10)
title('Population Map')
colormap gray

figure
pcolor(seafloor)
title('Population Map')
colormap pink
shading flat

figure
contour(topo,10)
title('Topographic Map')
colormap gray

figure
mesh((seafloor*.2)+10)
hold
pcolor(seafloor)
hidden off
colormap copper
shading flat

figure
surfl(seafloor,[-10 50])
colormap copper
shading flat
view(45,70)
axis([0 g 0 g 0 100])

%notification sound
Fs=10000;
x=sin(2*pi*4000/Fs*(0:1500));
sound(x,Fs)