function [x, fx, cnt] = NeldMead(fun,x,xtol,ftol,maxit)
prikaz = 1; pauza = 0;
disp('Prikazuje se svaki korak')
disp('za ukidanje prikaza postaviti: prikaz=0')
pause

n = prod(size(x));
maxit = maxit*n;

% Set up a simplex near the initial guess.
xin = x(:); % Force xin to be a column vector
v = xin;    % Place input guess in the simplex! (credit L.Pfeffer at Stanford)
fv = feval(fun,v);

% inicijalini simplex
for j = 1:n
   y = xin;
   if y(j) ~= 0
      y(j) = 1.2*y(j);
   else
      y(j) = 0.00025;
   end
   v = [v y];
   fv = [fv  feval(fun,y)];
end
[fv,j] = sort(fv);
v = v(:,j);
cnt = n+1;

if prikaz
   clc
   format compact
%   format short e
   home
   disp(' ')
   cnt
   disp('initial ')
   v
   fv
   if n==2, plot([v(1,:) v(1,1)],[v(2,:) v(2,1)]), drawnow, end
   if pauza, pause, end
end

alpha = 1;  beta = 1/2;  gamma = 2;
onesn = ones(1,n);
ot = 2:n+1;
on = 1:n;

% Iterate until the diameter of the simplex is less than tol.
while cnt < maxit
    if max(max(abs(v(:,ot)-v(:,onesn)))) <= xtol & ...
           max(abs(fv(1)-fv(ot))) <= ftol
        break
    end

    % One step of the Nelder-Mead simplex algorithm

    vbar = (sum(v(:,on)')/n)';
%    if prikaz, disp('Teziste'), disp(vbar), end
    vr = (1 + alpha)*vbar - alpha*v(:,n+1);
    fr = feval(fun,vr);  cnt = cnt + 1;
    vk = vr;  fk = fr; how = 'reflect ';
    if fr < fv(n)
        if fr < fv(1)
            ve = gamma*vr + (1-gamma)*vbar;
	    fe = feval(fun,ve);  cnt = cnt + 1;
            if fe < fv(1)
                vk = ve; fk = fe;
                how = 'expand  ';
            end
        end
    else
        vt = v(:,n+1); ft = fv(n+1);
        if fr < ft
            vt = vr; ft = fr;
        end
        vc = beta*vt + (1-beta)*vbar;
	fc = feval(fun,vc);  cnt = cnt + 1;
        if fc < fv(n)
            vk = vc; fk = fc;
            how = 'contract';
        else
            for j = 2:n
                v(:,j) = (v(:,1) + v(:,j))/2;
		fv(j) = feval(fun,v(:,j));
            end
            vk = (v(:,1) + v(:,n+1))/2;
   	    fk = feval(fun,vk);  cnt = cnt + n;
            how = 'shrink  ';
        end
    end
    v(:,n+1) = vk;
    fv(n+1) = fk;
    [fv,j] = sort(fv);
    v = v(:,j);

    if prikaz
        disp(' ')
        cnt
        disp(how)
        v
        size(v)
        fv
        if n==2, plot([v(1,:) v(1,1)],[v(2,:) v(2,1)]), drawnow, end

        %if pauza
          %keyboard
        %end
    end
end
x(:) = v(:,1);
if prikaz, format, end
fx=min(fv);
if cnt==maxit
    if prikaz
        disp(['Upozorenje: Maksimalan broj iteracija (', ...
               int2str(maxit),') je premasen']);
    end
end

*****************************
function demoNelderMead(fun,x0,xtol,ftol,maxit,x1span,x2span)
clc
optdemop(fun,x1span,x2span);
NeldMead(fun,x0,xtol,ftol,maxit)
***************************
function optdemop(f,x,y)
clc
disp(['Funkcija: ', f] )
hold off
clf
%
% racunanje vrednosti f-je 2 promenljive
%
[x,y] = meshgrid( x, y );
z=zeros(size(x));
pomoc = [x(:),y(:)]';
p = zeros(1,length(pomoc));
for i=1:length(pomoc), p(i)=feval(f,pomoc(:,i));, end
z(:) = p;
%
% 3-D
%
figure(1)
surf(x,y,z)
 shading interp;
colormap jet;
xlabel('x1'), ylabel('x2'), zlabel('f'), title('3D dijagram')
disp('Sa view(...) se moze menjati pogled na f-ju'),pause
%keyboard
%
% crtanje kontura funkcije
%
%clf
%contour(x,y,z,30,'k--')
%surf(x,y,z)

figure(2)
contour3(x,y,z,30)
surface(x,y,z,'EdgeColor','none','FaceColor','none')
%shading interp;
grid off
view(0,90)
 xlabel('x1'), ylabel('x2'), title('Konturni dijagram')
drawnow
hold on

dx=x(1,2)-x(1,1);
dy=y(2,1)-y(1,1);
[px,py] = gradient(z,dx,dy);
quiver(x,y,px,py)