function [f,x,nf,popout,popvelout] = PSO_Bound (fname,D,XVmin,XVmax,y,NP,itermax,omega,eta1,eta2,vmax,popin,popvelin)

if nargin == 4
    y=[];
    NP=20;
    itermax=1000;
    omega=0.65;
    eta1=2.0;
    eta2=2.0;
    vmax=.5;
end

%----pop is a matrix of size NPxD. It will be initialized-------------
%----with random values between the min and max values of the---------
%----parameters

if nargin < 12
    popin = zeros(NP,D); %initialize pop
    for i=1:NP
        popin(i,:) = XVmin + rand(1,D).*(XVmax - XVmin);
    end
end

%----popvel is a matrix of size NPxD. It will be initialized-------------
%----uniformly randomly within hyperparallelipiped-----------------------
%----of scale vmax of the search space when not passed as input
if nargin < 13
    popvelin = zeros(NP,D); %initialize popvel
    for i=1:NP
        popvelin(i,:) = (XVmin + rand(1,D).*(XVmax - XVmin))*vmax;
    end

end



personalbestVAL(1:NP)= Inf;
personalbestPOS = zeros(NP,D);
globalbestVAL = Inf;
val = zeros(1,NP);
nf=0;

%----pop is a matrix of size NPxD. It will be initialized-------------
%----with random values between the min and max values of the---------
%----parameters or passed as inputs to the optimiser---------

pop=popin;
popvel=popvelin;

%----Evaluate the objective function for each population member----%
for j=1:itermax
    %Prints Output every ten iterations
    if mod(j,itermax)==0
        fprintf(1,'Iteration: %d,  Best: %f,  om: %f,  et1: %f,  et2:%f, NP: %d\n',j,globalbestVAL,omega,eta1,eta2,NP);
        for n=1:D
            fprintf(1,'best(%d) = %f\n',n,globalbestPOS(n));
        end
    end
    %Updates the swarm velocities and positions
    for i=1:NP
        val(i) = feval(fname,pop(i,:),y);
        nf=nf+1;
        if val(i) < personalbestVAL(i)
            personalbestVAL(i)=val(i);
            personalbestPOS(i,:)=pop(i,:);
        end
        if val(i) < globalbestVAL
            globalbestVAL=val(i);
            globalbestPOS=pop(i,:);
        end
        r1=rand;
        r2=rand;
        popvel(i,:)=omega*popvel(i,:)+eta1*r1*(personalbestPOS(i,:)-pop(i,:))+eta2*r2*(globalbestPOS-pop(i,:));
        pop(i,:)=pop(i,:)+popvel(i,:);

        %Takes care that boundary constraints are enforced
        for i=1:NP
            for j=1:D
                if (pop(i,j) > XVmax(j)) | (pop(i,j) < XVmin(j))
                    pop(i,j) = XVmin(j) + rand * (XVmax(j) - XVmin(j));
                end
            end
        end
        x=globalbestPOS;
        f=globalbestVAL;
    end
end
    popout = personalbestPOS;
    popvelout = popvel;