% use graph-transform to compute manifold !

clear all;

a = 1.4;
b = 0.3;

phi_p = @(x)([ 1+x(2)-a*x(1).^2;b*x(1)]);

p(1) = (-(1-b)+sqrt( (1-b)^2+4*a))/(2*a);
p(2) = b*p(1);

% Compute in new coordinates !

phi =  @(x)( phi_p(x+p')-p' );


% Diagonalize

Dphi = Dapprox(phi,[0;0]);
[V,D]=eig(Dphi);

% < 1 => Stable manifold
% > 1 => Unstable manifold 
% Both manifolds are one dimensional

as = D(1,1);
au = D(2,2);

% Mapping to EU

u = 0:0.00001:0.001;
s = zeros(size(u));

% Construct new system in Es,Eu - coordinates !

W = inv(V);

psiS = @(s,u)(W(1,:)*phi( V*[s;u]));
psiU = @(s,u)(W(2,:)*phi( V*[s;u]));

us = u

for j=1:10

  disp('ITERATION');
  
  for i = 1:length(u)
    

    
    for k=1:100
      
      dus  = Dapprox( @(x)(psiU(sigma(s,u,x),x)-u(i)),us(i))\( ...
	  psiU(sigma(s,u,us(i)),us(i))-u(i));
      us(i) = us(i) -dus;
    
      norm(dus);
    
      if(norm(dus) < 1e-6)
	disp(strcat('Newton converged! :',num2str(i),'/',num2str(k)));
	
	s(i) = psiS( sigma(s,u,us(i) ),us(i) );
	
	break;
      end
    end
  end
end

for j=1:length(u)

  Y(j) = sigma(s,u,u(j));
  
end


Xu = V*[Y ;u];

%figure;
plot(Xu(1,:),Xu(2,:),'k*')
hold on;



%set(gca,'Xlim',[-0.2 0.2]);
%set(gca,'Ylim',[-0.2 0.2]);