function pts=Wu_g(p,nits,f)

% 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

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

A=Df(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

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)); 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));
      kp=floor(dz/dmax); 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 ,
      d=norm(z(:,i)-z(:,i+1))+norm(z(:,i+1)-z(:,i+2));
      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)); drawnow, nz
   y=z(:,1:nz); ny=nz; pts=[pts,z(:,2:nz)];
end; hold off
