function A4

x = linspace(0,1,21);               % Stützstellen
y = exp(x);                         % Funktionswerte
y_eps = y + randn(size(x))/1000;    % Gestörte Werte

xx = linspace(0,1,200);             % Auswertungsstellen
yy = exp(xx);

plot(xx,yy,'k');                    % Originale e-Funktion plotten
hold on;

plot(x,y,'.k');                     % Ungestörte Werte plotten
plot(x,y_eps,'.r');                 % Gestörte Werte plotten
legend('f(x)=e^x','y_i','y_i^\epsilon',2);
drawnow;
disp('Bitte warten...');

tic; flops(0);
P = Neville_Aitken(x,y,xx);         % Ungestörte Interpolation
fprintf('Zeit: %0.2fs    Flops: %i\n',toc, flops);
plot(xx,P,'g');
legend('f(x)=e^x','y_i','y_i^\epsilon','P mit P(x_i)=y_i',2);
drawnow;

tic; flops(0);
P_eps = Neville_Aitken(x,y_eps,xx); % Gestörte Interpolation
fprintf('Zeit: %0.2fs    Flops: %i\n',toc, flops);
plot(xx,P_eps,'r');
legend('f(x)=e^x','y_i','y_i^\epsilon','P mit P(x_i)=y_i','P mit P(x_i)=y_i^\epsilon',2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function f = Neville_Aitken(x,y,xj)
% Neville-Aitken-Algorithmus für die Stützstellen x
% und den Funktionswerten y. Dieser Algorithmus wertet
% das Interpolationspolynom durch die Punkte (x,y)
% an den Stellen xj aus.
%
% Eingabe: x  Stützstellen
%          y  Funktionswerte, die interpoliert werden
%          xj Stellen, an denen das Polynom ausgewertet wird
%
% Ausgabe: f  = P(xj), interpolierte Werte

n = length(x);
f = zeros(size(xj));

for j=1:length(xj)
   P = y;
   
   for k=2:n
      for i=n:-1:k
         P(i) = P(i) + (xj(j)-x(i)) / (x(i)-x(i-k+1)) * (P(i)-P(i-1));
      end;
   end;
   
   f(j) = P(n);
end;

