function [Tmax Vmax Ind Wmax] = getMinMax(T, Y, flag)

  % Get accurate voltage trace max or min using 4th order polynomial interpolation
  % Parameters are: T (time array) Y (phase space solution) flag (=1 for max, =-1 for min)
  %      Oh & Matveev, 2010

    global pf;
    ds = diff( sign( diff(Y(:,1)) ) );
    Ind = find(ds == -2*flag) + 1;    % Locations (indices) of critical points
    Tmax = zeros(size(Ind));
    Vmax = zeros(size(Ind));
    Wmax = zeros(size(Ind));

    for i = 1 : length(Ind)
        ind = Ind(i);
        pf = polyfit(T(ind-2:ind+2)-T(ind),Y(ind-2:ind+2,1)-Y(ind,1), 4);
        Tm =  fsolve( @(x) 4*pf(1)*x.^3 + 3*pf(2)*x.^2 + 2*pf(3)*x + pf(4), 0);  % Corrected time of voltage max/min
        Tmax(i) = T(ind) + Tm;
        Vmax(i) = Y(ind,1) + pf(1) * Tm^4 + pf(2)*Tm^3 + pf(3) * Tm^2 + pf(4)*Tm + pf(5);
        Wmax(i) = Y(ind,2);
    end
end
