function bonobeau(N,X,alphaa,gammaa,gen,sd,css,sync,moviee)
% Generates Bonabeau swarm.
%
% PARAMETERS
%   N       =   Number of agents
%   X       =   Size of Playground
%   alphaa  =   Swarm Parameter = 1 - betaa
%   gen     =   Generations (number of iterations)
%   nowall  =   If 1, the Playground is modulo.  
%               If 0, there is reflection off of walls
%               Else, the agent disappears
%   sd    =   If 0, the initilization is random
%               Else, a centered Gaussian with standard deviation 'sd'
%   css   =   If 1, there is a constant step size of gammaa.
%               Else, step size is in accordance to dynamic equation
%   sync  =    If 1, the update is synchronous
%              If 0, the agents are shuffled each epoch
%               Else, agents are chosen randomly
%   moviee =   If 1, make an avi movie.  Otherwise, no.
%
% BLACK AND RED TRIOS ARE COMMENTED OUT, BUT CAN BE INSERTED
% Two agent trios are followed, one black and one red.
%       The circle 'o' follows the 'x' and the '+'
%       For alphaa > 0.5, the 'o' is trying to get closer to the '+' than
%       the 'x'
%
%                           
close all;
decayy=0.995;
if moviee==1
    fig=figure;
    mov = avifile('example.avi')
end    
Playground=zeros(X,X);
Playground(1,1)=1; % Used for color reasons.
if sd==0 
    p=(X+1)*rand(N,2);
else
    p=X/2*ones(N,2)+sd*randn(N,2);
    %p=max(zeros(size(p)),p);
    %p=min(X*ones(size(p)),p);
end
betaa=1-alphaa;
ii=zeros(N,2);
pdot=zeros(N,2);
for n=1:N
    for m=1:N
        ii(n,1)=ceil(N*rand);
        if ii(n,1)~=n; break; 
        end
    end
    for m=1:N
        ii(n,2)=ceil(N*rand);
        if and(ii(n,2)~=n, ii(n,2)~=ii(n,1)); break; 
        end    
    end
end
for g=1:gen
    if sync==0
        [Y,I] = sort(rand(N,1));
    end
    for m=1:N
        if sync==1
            n=m;
        elseif sync==0;
            n=I(m);
        else
            n=ceil(N*rand);
        end
        pdot(n,:)=-gammaa*(p(n,:)-(alphaa*p(ii(n,1),:)+betaa*p(ii(n,2),:)));
        if css==1
            denn=pdot(n,1)^2+pdot(n,2)^2;
            p(n,:)=p(n,:)+gammaa*pdot(n,:)/sqrt(denn);
        else
            p(n,:)=p(n,:)+pdot(n,:);
        end
    end
    Playground=decayy*Playground;
    imagesc(Playground);hold on;
    for n=1:N 
        plot(p(n,1),p(n,2),'w.');hold on
        plot(p(n,1),p(n,2),'ko');hold on
        yyy=round(p(n,2));xxx=round(p(n,1));
        bitt=or(or(xxx<1,yyy<1),or(yyy>X,xxx>X));
        if bitt==0  
            Playground(yyy,xxx)=1;
        end
    end
    title(g);
    F = getframe;
    if moviee==1
        mov = addframe(mov,F);
    end
    hold off
end  % gen
if moviee==1
    mov = close(mov);
end

Resources

e website templates.