The Invasive Weed Optimization Algorithm: An Implemenation

I am working on implementation of the invasive weed optimization (IWO) algorithm since one year or so. The selection of this algorithm has linkage to my family business background which is very closely related with manufacture and maintain tools to support agriculture, vegetation. The original IWO algorithm was proposed in a paper by Mehrabian and Lucas (2006) with title, "a novel numerical optimization algorithm inspired from weed colonization." I have made a number of attempts to write a program code in MATALB. The basic working version is recently completed and pasted below. The algorithm has four major steps: (1) initialization, (2) reproduction, (3) spatial dispersal, and (4) competitive exclusion. The program code is pasted below.

% created by dr tariq javid as of 02-december-14
% last updated on 21-november-15
% f1 = sphere, f2 = griewant, f3 = rastrigin
clear all, clc
% --- user input
init_pop = input('number of initial population [10]: ');
iter_max = input('maximum number of iterations [100]: ');
dimn_prb = input('problem dimension [2]: ');
plnt_max = input('maximum number of plant population [15]: ');
seed_max = input('maximum number of seeds [5]: ');
seed_min = input('minimum number of seeds [0]: ');
nonl_mid = input('nonlinear modulation index [3]: ');
stdv_ini = input('initial value of standard deviation [3]: ');
stdv_fin = input('final value of standard devition [0.001]: ');
srch_ini = input('a. initial search area - lower [-40]: ');
srch_fin = input('b. initial search area - upper [-30]: ');
% --- default value
if isempty(init_pop), init_pop = 10   ; end
if isempty(iter_max), iter_max = 100  ; end
if isempty(dimn_prb), dimn_prb = 2    ; end
if isempty(plnt_max), plnt_max = 15   ; end
if isempty(seed_max), seed_max = 5    ; end
if isempty(seed_min), seed_min = 0    ; end
if isempty(nonl_mid), nonl_mid = 3    ; end
if isempty(stdv_ini), stdv_ini = 3    ; end
if isempty(stdv_fin), stdv_fin = 0.001; end
if isempty(srch_ini), srch_ini = -40  ; end
if isempty(srch_fin), srch_fin = -30  ; end
% --- function definition
f1 = @(x,y) x.^2 + y.^2; % for default values
% f2 = @(x,y) (x.^2 + y.^2)/4000 - cos(x/sqrt(2)) .* cos(y/sqrt(2)) + 1;
% f3 = @(x,y) 20 + x.^2 + y.^2 - 10 * (cos(2*pi*x) + cos(2*pi*y));
% --- initialize weeds as random numbers
weed_ini = srch_ini + (srch_fin - srch_ini).*rand(init_pop,dimn_prb)
% --- preparation of next step
temp_one = (stdv_ini - stdv_fin); temp_two = iter_max .^ nonl_mid ;
weed_pop = weed_ini;
figure, axis([-50 10 -50 10]), hold on
plot(weed_ini(:,1),weed_ini(:,2),'r^')
for iter_idx = 1:iter_max % main for loop
    temp_thr = ((iter_max - iter_idx) .^ (nonl_mid)) / temp_two;
    iter_now = temp_thr .* temp_one + stdv_fin;
    % --- evaluate fitness of weeds
    for weed_idx = 1:length(weed_pop)
        comp_fit(weed_idx) = f1(weed_pop(weed_idx,1),weed_pop(weed_idx,2));
    end
    % --- determine number of seeds as per fitness
    fit_max = min(comp_fit); fit_min = max(comp_fit);
    slope = (seed_max - seed_min) / (fit_max - fit_min);
    seed_cnt = zeros(size(comp_fit));
    for seed_idx = 1:length(weed_pop)
        seed_tmp = floor(slope * (comp_fit(seed_idx) - fit_min) + seed_min);
        seed_cnt(seed_idx) = seed_tmp;
    end
    % --- spatial dispersal
    for weed_idx = 1:length(weed_pop)
        if seed_cnt(weed_idx) ~= 0 % to avoid nan and inf unexpected termination
            weed_tmp = iter_now .* randn(seed_cnt(weed_idx),dimn_prb);
            for seed_idx = 1:seed_cnt(weed_idx)
                weed_tmp(seed_idx,1) = weed_tmp(seed_idx,1) + weed_pop(weed_idx,1);
                weed_tmp(seed_idx,2) = weed_tmp(seed_idx,2) + weed_pop(weed_idx,2);
            end
            weed_pop = cat(1,weed_pop,weed_tmp);
        end
    end
    % --- evaluate fitness of new population
    for weed_idx = 1:length(weed_pop)
        comp_fit(weed_idx) = f1(weed_pop(weed_idx,1),weed_pop(weed_idx,2)); 
    end
    % --- competitive exclusion
    while length(weed_pop) > plnt_max
    % --- i is index of maximum value i.e. worse fit
       [~, i] = max(comp_fit);
       comp_fit(i)   = [];    % remove maximum value from fitness
       weed_pop(i,:) = [];    % remove worse fit weed from population
    end
    % --- plot result for each iteration
    plot(weed_pop(:,1),weed_pop(:,2),'b*');
end
weed_pop
for weed_idx = 1:length(weed_pop)
    comp_fit(weed_idx) = f1(weed_pop(weed_idx,1),weed_pop(weed_idx,2));
end
[~, i] = min(comp_fit);
weed_idx = i
weed_sol = [weed_pop(i,1),weed_pop(i,2)]
weed_fit = comp_fit(i)
plot(weed_pop(i,1),weed_pop(i,2),'r^');
legend('initial,final','progress',2)
grid on, hold off