You are an expert in writing code in many different languages. Your goal is to perform code completion for the following code, keeping in mind the rest of the code and the file meta data.
    Metadata details:
    filename: max.m
    path of file: ./repos/mtex/SO3Fun/@SO3Fun
    
    the code of the file until where you have to start completion: function [values,modes] = max(SO3F,varargin)
% global, local and pointwise maxima of functions on SO(3)
%
% Syntax
%   [v,pos] = max(SO3F) % the position where the maximum is atained
%
%   [v,pos] = max(SO3F,'numLocal',5) % the 5 largest local maxima
%
%   SO3F = max(SO3F, c) % maximum of a rotational functions and a constant
%   SO3F = max(SO3F1, SO3F2) % maximum of two rotational functions
%   SO3F = max(SO3F1, SO3F2, 'bandwidth', bw) % specify the new bandwidth
%
%   % compute the maximum of a multivariate function along dim
%   SO3F = max(SO3Fmulti,[],dim)
%
% Input
%  SO3F, SO3F1, SO3F2 - @SO3Fun
%  SO3Fmulti          - a multivariate @SO3Fun
%  c                  - double
    
function on grid
  f = eval(SO3F,S3G,varargin{:});

  % take only local minima as starting points
  ind = isLocalMinD(-f,S3G,1.5*res0,varargin{:});
  modes = S3G(ind);
  values = f(ind);

  % the total walking distance
  sumOmega = zeros(size(modes));

  numLocal = get_option(varargin, 'numLocal',1);

  % local refinement
  res = res0;
  while res > targetRes
    res = res / 1.25;

    % neighborhood search
    S3Glocal = localOrientationGrid(crystalSymmetry,crystalSymmetry,2*res,'resolution',res/2);
    newModes = (S3Glocal * modes).';
    f = eval(SO3F,newModes,varargin{:});
  
    if numLocal == 1

      [values,id] = max(f(:));
      modes = newModes(id);
        
    else
      [values,id] = max(f,[],2);
      modes = newModes(sub2ind(size(newModes),(1:length(modes)).',id));

      % keep track of the walking distance
      omega = S3Glocal.angle;
      sumOmega = sumOmega + omega(id);

      % maybe we can reduce the number of points a bit
      [~,~,I] = unique(modes, 'tolerance', 1.5*res0,'noSymmetry');
      modes = normalize(accumarray(I,modes));
      values = accumarray(I,values,[],@mean);
      sumOmega = accumarray(I,sumOmega,[],@min);

      % consider only points that did not walked too far
      values(sumOmega>5*res0) = [];
      modes(sumOmega>5*res0) = [];
      sumOmega(sumOmega>5*res0) = [];
    end