% tht1 !


% define functions needed !
f = @(x)( x(1)^2 + x(1) + x(2)^2 - x(2) - 1 );
Df = @(x)( [2*x(1)+1 , 2*x(2)-1] );


% plot stuff
r = sqrt(3/2);
t = 0:0.01:2*pi;
y = sin(t)*r+1/2;
x = cos(t)*r-1/2;

plot(x,y,'r-');
hold on;

axis square;



while(1)

  % wait for input
  
  [x,y,BUTTON] = ginput(1);
  
  if(BUTTON == 3)
    break;
  end
  
  
  x0 = [ x;y]
  
  % solve !
  [useless,X] = gnewton(x0,f,Df);
  [useless,X2] = gnewton(x0,f,Df,@luifun);

  % plot iterates 
  
  plot(X(1,:),X(2,:),'r*-');
  
  abu = -(2*X(1,:)+1)./(2*X(2,:)-1);
  abu(find(isinf(abu))) = 1;
  quiver(X(1,:),X(2,:),X(1,:),abu.*X(1,:));
  
  
  plot(X2(1,:),X2(2,:),'ko--');
  
end
