function [] = prog1main()
mode(0);
ieee(1);
  nn=[5,10,20];
  for k=1:3
  n=nn(k);
  for i=1:n+1
    x(i)=-5+(i-1)*10/n;         //stuetzstelle
  end
  for i=1:n
    y(i)=-5+5/n+(i-1)*10/n;     //Mittelpunkte y(i)
  end

// bildX() bildet die Matrix x der Newtonsche Interpolation.
// urf(vektor) loest den exakte Wert der f(x) in der Stelle x.
// rA(matrix A,vektor B) loest den Wert x
  a=rA(bildX(x),urf(x));
//p(vektor a, vektor y) loest die numerische Loesung bei Newtonsche Interpolation
  pp=p(a,y);
  ff=urf(y);
  verg=0;
  for i=1:length(pp) //berechnet den maximalen Fehler
      if (verg<abs(pp(i)-ff(i)))
         verg=abs(pp(i)-ff(i));
      end
  end
  disp('Newtonsche Interpolation')
  disp(p(a,y))  
  disp(urf(y))
  disp('maximalen Fehler:')
  disp(verg)
  
  
  disp('natuerliche kubisch Spline')
  //bildB() rechnet den Vektor b ,was danach braucht.
  b=bildB(urf(x));
  //bildXX() rechnet die Matrix A, was danach braucht
  A=bildXX(n);
  //rS(matrix,vektor) loest das Gleichungssystem und bekommt man den Vektor s,was Spline interp braucht.
  S=rS(A,b);
  s(1)=S(1);
  s(2)=S(n+1);
  //p2(vektor s, vektor y) loest die numerische Loesung bei der natuerliche kubische Spline
  pp=p2(s,y);
  verg=0;
  for i=1:length(pp)  //berechnet den maximalen Fehler 
      if (verg<abs(pp(i)-ff(i)))
         verg=abs(pp(i)-ff(i));
      end
  end
  disp(p2(a,y))  
  disp(urf(y))
  disp('maximalen Fehler:')
  disp(verg)  
 end
 
endfunction

function f=urf(x)
   n=length(x);
   for i=1:n
     f(i)=1/(x(i)^2+1);
   end
endfunction

function p=p(a,xx)
 m=length(xx);
 n=length(a);
 for k=1:m
   for i=1:n
      x(i)=-5+(i-1)*10/(n-1);
   end
   p(k)=0;
   for i=1:n
      mal=1;
      for j=1:(i-1)
        mal=mal*(xx(k)-x(j));
      end
      p(k)=p(k)+a(i)*mal; 
   end 
 end
endfunction

function pp=p2(s,x)
     m=length(x);
     for i=1:m
         pp(i)=urf(-5)+s(1)*(x(i)-(-5))+(3*(urf(5)-urf(-5))/10^2-(s(2)+2*s(1))/10)*(x(i)-(-5))^2+((-2)*(urf(5)-urf(-5))/10^3+(s(2)+s(1))/10^2)*(x(i)-(-5))^3
     end
endfunction     
     
function X=bildX(x)
    n=length(x);
    for i=1:n
      for j=1:i
          mal=1;
          for k=1:j-1
            mal=mal*(x(i)-x(k));
          end
          X(i,j)=mal;
      end
    end
endfunction

function X=bildXX(n)
         h=10/(n);
         X(1,1)=2;
         X(1,2)=1;
         for i=2:n
             X(i,i)=4*h;
             X(i,i-1)=h;
             X(i,i+1)=h;
         end
         X(n+1,n)=1;
         X(n+1,n+1)=2;   
endfunction

function b=bildB(f)
       n=length(f);
       h=10/(n-1);
       b(1)=3*(f(2)-f(1))/h;
       for i=2:n-1
          b(i)=3*(f(i+1)-f(i-1));
       end
       b(n)=3*(f(n)-f(n-1))/h;
endfunction

function a=rA(x,b)
   n=length(b);
   for i=1:n 
      summ=0;
      for j=1:i-1
          summ=summ+x(i,j)*a(j);
      end
      a(i)=(b(i)-summ)/x(i,i);
   end 
endfunction

function x=rS(A,b)
  n=length(b);
  x(1:n)=0;
  for k=1:300
    N(1:n)=0;
    M(1:n)=0;
    for i=1:n
        N(i)=(1/A(i,i))*b(i);
        for j=1:i-1
            M(i)=M(i)+(-A(i,j))*x(j);
        end
        for j=i+1:n
            M(i)=M(i)+(-A(i,j))*x(j);
        end
        M(i)=M(i)*(1/A(i,i));
        x(i)=M(i)+N(i); 
     end
  end 

endfunction
