
function pts = Wu_g2(p,nits,f,c)

% Input:
%  p        fixed point (a column vector)
%  nits     number of required mappings of initial segment
%  f        string that names the function whose unstable 
%           manifold of the fixed point is sought

if(nargin < 4)
  c = 'b-';
end


dmax=0.01; dim=length(p);
zers=zeros(dim,1000); z=zers; zmax=1000;

A=Dapprox(f,p); [V,D]=eig(A);
[lambda,j]=max(abs(diag(D))); if D(j,j)<0 , return; end
v=V(:,j)/norm(V(:,j)); % this gets one branch of the manifold
                       % for the other, change the sign

y(:,1)=p+dmax*v/lambda; y(:,2)=feval(f, y(:,1)); ny=2;
pts=[p,y(:,1:2)]; np=3; 
plot(y(1,1:ny),y(2,1:ny),c); drawnow, hold on

for k=1:nits ,
   z(:,1)=y(:,ny); nz=1;
   for i=2:ny , 
      v=feval(f,y(:,i)); dz=norm(v-z(:,nz));
      kcurv=0; if nz>1; d1=v-z(:,nz); d2=z(:,nz)-z(:,nz-1);
        kcurv=40*round(1-d1'*d2/(norm(d1)*norm(d2))); end 
      kp=max(ceil(dz/dmax),kcurv); nz=nz+kp+1; dy=y(:,i)-y(:,i-1);
           if nz>zmax, zmax=zmax+1000; z=[z,zers]; end,
      for j=1:kp , z(:,nz-j)=feval(f,y(:,i)-j/(kp+1)*dy); end
      z(:,nz)=v; end;

   i=1; while i<nz-1 ,
      d1=z(:,i)-z(:,i+1); d2=z(:,i+1)-z(:,i+2);
      n1=norm(d1); n2=norm(d2);
      d=max(n1+n2,0.1*(1-d1'*d2/(n1*n2)));
      if d<dmax , z(:,i+1:nz-1)=z(:,i+2:nz); nz=nz-1; 
      else i=i+1; end; end;
   plot(z(1,1:nz),z(2,1:nz),c); drawnow, nz
   y=z(:,1:nz); ny=nz; pts=[pts,z(:,2:nz)]; 
end; hold off