function [p,x,E1all,E2all,xall]=remez(f,a,b,x,maxiter,eps)

% Remez-Algorithmus zur Polynom-Approximation
% der Funktion f:[a,b]->R mit Polynomen vom
% Grad n
%
% Input:
% x:  Vektor, enthält die Startreferenz
%     n=length(x)-2 ist der Polynomgrad
% maxiter: Max. Zahl der Iterationen
% eps: Abbruch des Verfahrens bei Fehlergenauigkeit eps
%
% Output:
% p: Polynom der Bestapproximation, Ausgabe erfolgt als
%    Koeffizientenvektor p=[c0,c1,c2,...,cn] mit den
%    Koeffizienten der Newton-Polynome  1, x-x(1), ...
% x: Vektor der letzten Referenz (Approx. einer Extremal-
%      alternante)
% E1: obere Schranke des Approx.grades für alle Iterierten
% E2: untere Schranke des Approx.grades für alle Iterierten
% xall: Matrix, enthält die Referenzen der gesamten Iteration als Zeilen
% 
% Bei erfolgreichem Ende ist E1-E2<eps

x=x(:)'; % stelle einen Zeilenvektor her
x=sort(x);
xall=x;
n=length(x)-2; % Polynomgrad
tn=10001; % Diskretisierung für Berechnung des Maximums
t=linspace(a,b,tn); % Koordinaten zur Berechnung 
     % von ||f-p||

y1=(-1).^[0:n+1]; % Vektor der Vorzeichen +-1
k=0; % Iterationszähler
EE=eps+1;
E1all=[];
E2all=[];

while ((EE>eps) & (k<maxiter))
    % löse diskretes Approx.problem zur Referenz x
    [c1,w1]=divdiff(x,y1);
    [c2,w2]=divdiff(x,feval(f,x));
    d=c2/c1;
    p=w2-d*w1; 
    p=p(1:end-1); % diskrete Best-Approx. in Newton-Form 
    E2=abs(d); % Fehler der diskreten Approx.
    % hier beginnt der Austauschschritt
    z1=feval(f,[t,x]);
    z2=newton_eval(p,x(1:end-1),[t,x]);
    [E1,i]=max(abs(z1-z2)); % max. Fehler bei t(i(1))
    xneu=t(i(1));
    EE=E1-E2;
    % Plot der diskreten Best-Approx.
%    figure
    clf
    plot(t,z1(1:tn)) % plot f
    hold on
    plot(t,z2(1:tn),'r') % plot p
    x_bars=[x;x;NaN*ones(1,n+2)];
    x_bars=x_bars(:)';
    y_bars=[z1(tn+1:end);z2(tn+1:end);NaN*ones(1,n+2)];
    y_bars=y_bars(:)';
    plot(x_bars,y_bars,'r') % markiere Referenz  
    plot([xneu,xneu],[z1(i(1)),z2(i(1))],'mv')
    % Austausch
    if ((EE>eps) & (k<(maxiter-1)))
        % finde auszutauschendes x(j)
        j=find((xneu-x)>0);
        if length(j)>0
            jj=j(end);
            xj=x(j(end));
            z1=feval(f,[xj,xneu]);
            z2=newton_eval(p,x(1:end-1),[xj,xneu]);
            k=((z1(1)-z2(1))*(z1(2)-z2(2)))<0;
            x(jj+k)=xneu;
        else
            % xneu ist kleiner als das alte x(1)
            x(1)=xneu;
        end
        xall=[xall;x];
    end
    E1all=[E1all;E1];
    E2all=[E2all;E2];    
    k=k+1;
    pause
end