% Das Programm INTERPOL interpoliert die Funktionen
% 
%   f(t) = 1/(1+25t^2)   und   g(t) = sqrt(|t|)
%
% im Intervall [-1,1] an äquidistanten bzw. 
% Tschebyscheff-Stützstellen x(i), i=0,...,n
% für n=2,4,6,...,20.

film=0;
k=menu('Interpolation mittels Newtonschema',...
       'Äquidistante Stützstellen',...
       'Chebyshov-Stützstellen',...
       'Ende');
while k<3
      clear y x;
      v=2:2:20;     
      lf=menu('Zu interpolierende Funktion',...
	      '1/(1+25x^2)',...
	      'sqrt(|x|)');
      for n=v
	  if k==1
	     t =-1:(2/n):1;
	     if lf~=1	     
	     f=sqrt(abs(t));
	     else
	     f=1./(1+25*t.^2);
	     end
	     x = t(1):((t(2)-t(1))/21.0):t(2);
	     for i = 2:length(t)-1
	         x = [x(1:length(x)-1),t(i):((t(i+1)-t(i))/21.0):t(i+1)];
	     end % for i
	     for ii=1:length(t)
             dd(ii)=f(ii);
                for kk=1:ii-1
                   dd(ii) = (dd(ii) - dd(kk))/(t(ii)-t(kk));
                end % for kk
             end % for ii
             for ii=1:length(x)
		 mm=length(t);
                 y(ii)=dd(mm)*ones(1,length(x(ii)));
                 for jj=mm-1:-1:1
                    y(ii) = y(ii) .* (x(ii) - t(jj)) + dd(jj);
                 end % for jj
	     end  % for ii
	     if lf~=1	     
	     z=sqrt(abs(x));
	     else
	     z=1./(1+25*x.^2);
	     end	     
	     figure(1);
	     plot(x,z,'-g',x,y,':r');
	     w=axis;
	     axis([w(1:2) -5 5]);
	     title('Interpolationspolynom Teil i) - Stützstellen in blau');
	     text(0.5,4,['n=' num2str(n)]);
	     xlabel('t');
	     ylabel('y');
	     hold on;
	     plot(t,f,'o');
	     hold off;
	     if film 
	        MV=getframe;
		MV1(1:length(MV),j)=MV;
	     end;
	  end;	     
	  if k==2 
	     t = cos((2*(0:n)+1) * (pi/(2*(n+1))));
	     if lf~=1	     
	     f=sqrt(abs(t));
	     else
	     f=1./(1+25*t.^2);
	     end
	     x = t(1):((t(2)-t(1))/21.0):t(2);
	     for i = 2:length(t)-1
	         x = [x(1:length(x)-1),t(i):((t(i+1)-t(i))/21.0):t(i+1)];
	     end % for
	     for ii=1:length(t)
             dd(ii)=f(ii);
                for kk=1:ii-1
                   dd(ii) = (dd(ii) - dd(kk))/(t(ii)-t(kk));
                end % for kk
             end % for ii
	     for ii=1:length(x)
		 mm=length(t);
                 y(ii)=dd(mm)*ones(1,length(x(ii)));
                 for jj=mm-1:-1:1
                    y(ii) = y(ii) .* (x(ii) - t(jj)) + dd(jj);
                 end % for jj
	     end % for ii
	     if lf~=1	     
	     z=sqrt(abs(x));
	     else
	     z=1./(1+25*x.^2);
	     end
	     figure(2);
	     plot(x,z,'-g',x,y,'--r');
	     w=axis;
	     axis([w(1:2) -0.5 1.5]);
	     title('Interpolationspolynom Teil ii) - Stützstellen in blau');
	     text(0.5,1,['n=' num2str(n)]);
	     xlabel('t');
	     ylabel('y');
	     hold on;
	     plot(t,f,'o');
	     hold off;
	     if film
	        MV=getframe;
		MV2(1:length(MV),j)=MV;
	     end;
	  end % if
	  j=j+1;
      end % for
      if film
	 if k==1 
	    figure(1);
	    movie(MV1,1,1);
	 end % if
	 if k==2 
	    figure(2);
	    movie(MV2,1,1);
	 end % if
      end
      k=menu('Interpolation mittels Newtonschema',...
	     'Äquidistante Stützstellen',...
	     'Chebyshov-Stützstellen',...
	     'Ende');
end % while
