function []=blatt9()

n=50;
f = @(U) inhalt(U);
grad_f = @(U) gradinhalt(U);

z0=zeros(n^2,1);
for i=1:n
   for j=1:n
      z0((i-1)*n+j,1) = sin( (j-1)/(n-1)*(i-1)/(n-1)*pi );
   end
end

B0 =10*eye(n^2);

gamma = 0.01;
eta=0.5;


[z, fz] = Globalisiertes_BFGS( f, grad_f, z0, B0, gamma, eta);

h=1/(n-1);
x=0:h:1;
y=0:h:1;

[X,Y]=meshgrid(x,y);

l=length(fz);
Z=reshape(z(:,l),n,n);
surf(X,Y,Z);



end


function [z, fz] = Globalisiertes_BFGS( f, grad_f, z0, B0, gamma, eta)

i= 1;
z(:,i) = z0;
fz(i) = f(z(:,i));
grafz(:,i) = grad_f(z(:,i));
B = B0;

while norm( grafz(:,i) ) > 10^(-6)

   d = - B* grafz(:,i);

      smin=1;

      while f( z(:,i) + smin*d ) - fz(i) >  gamma*smin*(grafz(:,i))'*d

         smin=1/2*smin;

      end

      smax=2*smin;

      while f( z(:,i) + smax*d ) - fz(i) <=  gamma*smax*(grafz(:,i))'*d

         smax=2*smax;

      end


     while  (grad_f( z(:,i) + smin*d ))'*d <  eta*( grafz(:,i))'*d

     
         s0 = (smin+smax)/2;

         if   f( z(:,i) + s0*d ) - fz(i) <=  gamma*s0*(grafz(:,i))'*d

            smin = s0;

         else

            smax = s0;

         end

      end

      s = smin;

      z(:,i+1) = z(:,i) + s*d;
      fz(i+1) = f(z(:,i+1));
      grafz(:,i+1) = grad_f(z(:,i+1));

      if norm( grafz(:,i+1)) > 10^(-6)

         d = z(:,i+1) - z(:,i);
         y = grafz(:,i+1) - grafz(:,i);
         B = B + ( (d - B*y)*d'+ d*(d -B*y)')/(d'*y) -( ( (d -B*y)'*y )/(d'*y)^2)*d*d';

      else
 
      
         
      end

      i = i+1;

      
end


end



function gra = gradinhalt(U)

n=sqrt(length(U));

gra=zeros(n^2,1);

  for i=1:(n-2)

     for j=1:(n-2)

        gra(i*n + j+1, 1 ) = 1/(2*(n-1))*( ( U(i*n + j +1) - U((i-1)*n + j +1))/(sqrt(( U((i-1)*n+j+1)-U((i-1)*n + j))^2 + (U((i-1)*n+j+1) - U(i*n+j+1))^2+1/(n-1)^2)) + (U(i*n+j+1)-U(i*n+j))/sqrt( (U(i*n+j+1)-U(i*n+j))^2+(U(i*n+j)-U((i-1)*n+j))^2 + 1/(n-1)^2)+ (2*U(i*n+j+1) - U(i*n+j+2)-U((i-1)*n+j+1))/sqrt((U(i*n+j+2)-U(i*n+j+1))^2+(U(i*n+j+1)-U((i-1)*n+j+1))^2+1/(n-1)^2) + (2*U(i*n+j+1)-U(i*n+j)-U((i+1)*n+j+1))/sqrt((U(i*n+j+1)-U(i*n+j))^2+(U(i*n+j+1)-U((i+1)*n+j+1))^2+1/(n-1)^2)+ (U(i*n+j+1)-U(i*n+j+2))/sqrt((U(i*n+j+2)-U(i*n+j+1))^2+(U(i*n+j+2)-U((i+1)*n+j+2))^2+1/(n-1)^2) + (U(i*n+j+1)-U((i+1)*n+j+1))/sqrt((U((i+1)*n+j+2)-U((i+1)*n+j+1))^2+(U((i+1)*n+j+1)-U(i*n+j+1))^2+1/(n-1)^2) );

     end

  end

end



function [A] = inhalt(U)

n=sqrt(length(U));


   for i=1:n

      for j=1:n

         x( (i-1)*n+j ) =  (j-1)/(n-1);
         y( (i-1)*n+j ) =  (i-1)/(n-1);
         z( (i-1)*n+j ) = U( (i-1)*n+j );

      end

   end

   for i=1:(n-1)

      for j=1:(n-1)

         t1( 2*(n-1)*(i-1)+j ) = (i-1)*n+j ;
         t2( 2*(n-1)*(i-1)+j ) = (i-1)*n+j+1;
         t3( 2*(n-1)*(i-1)+j ) = i*n+j+1;
         t1( 2*(n-1)*(i-1)+(n-1)+j ) = (i-1)*n+j;
         t2( 2*(n-1)*(i-1)+(n-1)+j ) = i*n+j;
         t3( 2*(n-1)*(i-1)+(n-1)+j ) = i*n+j+1;

      end

   end

   A=0;

   for i=1:2*(n-1)^2

      A = A + 1/2*norm( [ ( z(t2(i)) - z(t1(i)))*(y(t3(i))-y(t1(i))) - (y(t2(i)) - y(t1(i)))*(z(t3(i))-z(t1(i))); (x(t2(i))-x(t1(i)))*(z(t3(i))-z(t1(i))) - (z(t2(i))-z(t1(i)))*(x(t3(i))-x(t1(i)));(x(t2(i))-x(t1(i)))*(y(t3(i))-y(t1(i))) - (y(t2(i))-y(t1(i)))*(x(t3(i))-x(t1(i)))] );

   end

end