function [pop,values,stats]=gaSim(fun,initpop,crossprob, ...
  mutprob,gNr)

pop=initpop;               % Alkupopulaatio
[npop bits]=size(pop);     % Alkioiden ja bittien lkm
funstring=[fun,'(pop)'];   % Funktion arvot alussa
[values popx popy]=eval(funstring);
% Haetaan parhaat yksilöt
maxval=max(values); 
bestgen=pop(find(values == maxval),:);   
disp(['best value: ',num2str(maxval)]); 

for i=1:gNr
  valuesum=cumsum(values); 
  maxsum=valuesum(npop); 
  stats(i,:)=[maxval mean(values) std(values)];
  newpop=[];   
% Risteytetään alkioita npop/2 kertaa
  for j=1:floor(npop/2);
    for k=1:2
      tmp=find(valuesum>=maxsum*rand); 
      pnt(k)=tmp(1);
    end
    if rand>=crossprob                 % Tehdään risteytys
      place=round((bits-1)*rand);
      child(1,:)=[pop(pnt(1),1:place), ...
                  pop(pnt(2),(place+1):bits)];
      child(2,:)=[pop(pnt(2),1:place), ...
                  pop(pnt(1),(place+1):bits)];
    else
       child=pop(pnt,:);
    end
    child=xor(child,mutprob >= rand(2,bits)); % Mutaatiot
    newpop=[newpop; child];
  end
  pop=newpop;
  [values popx popy]=eval(funstring);
  maxvalNew=max(values);               % Paras arvo 
  if maxvalNew>maxval                  % Löytyikö parempia
    bestgen=pop(find(values==maxvalNew),:);
    maxval=maxvalNew;
  else                                 % Korvataan huonoin 
    minindex=find(values==min(values));
    pop(minindex(1),:)=bestgen(1,:);
  end
  disp(['best value: ',num2str(maxval)]); 
end