h=plot([-pi,0,pi],0,'k.'); set(h,'markersize',40); hold on; 
	pdl = inline('[x(2);-9.81*sin(x(1))]','t','x');
	opt=odeset('RelTol',1e-5);
	for x1=[-9 -8 -7 1 2 3 4 5 6 7 8 9],
	    [t,x]=ode45(pdl,[0 5],[0;x1],opt);
	    h=plot(x(:,1),x(:,2),'k-');  set(h,'linewidth',2);
	    [t,x]=ode45(pdl,[0 -5],[0;x1],opt);
	    h=plot(x(:,1),x(:,2),'k-');  set(h,'linewidth',2);
	end
	axis([-3.2 3.2 -8 8])
	
	[t,x]=ode45(pdl,[0 -4],[3.14;0],opt);
	h=plot(x(:,1),x(:,2),'-','Color',[0.6 0.6 0.6]); set(h,'linewidth',4);
	
	[t,x]=ode45(pdl,[0 4],[3.14;0],opt);
	h=plot(x(:,1),x(:,2),'-','Color',[0.6 0.6 0.6]); set(h,'linewidth',4);
					

	n = 30;	tau = 0.1;	x = zeros(2,n);	x(:,1) = [0;4];
	for i=1:n,
	    x(:,i+1) = x(:,i) + tau*feval(pdl,0,x(:,i));
	end
	h=plot(x(1,:),x(2,:),'k.-'); set(h,'linewidth',2,'markersize',20);
					

	x0 = [0;4];	n = 30;	tau = 0.1;	x = zeros(2,n);	x(:,1) = x0;
	TOL = 1e-6;	N_max = 20;	d = size(x0,1);	x(:,1) = x0;
	for k = 1:n
	    x(:,k+1) = x(:,k);
	    for j = 1:N_max
	        % Newton step 
	        F = x(:,k+1) - (x(:,k) + tau * feval(pdl,0,x(:,k+1)));
	        DF = eye(d) - tau * feval(dpdl,0,x(:,k+1));
	        dx = DF\F;
	        x(:,k+1) = x(:,k+1) - dx;
	        if (max(abs(dx)) < TOL) break; end;
	    end;
	end;
	h=plot(x(1,:),x(2,:),'k.-');
	set(h,'linewidth',2,'markersize',20);
					

	x0 = [0;4]; n = 30; tau = 0.1; x = zeros(2,n);
	x(:,1) = x0; TOL = 1e-6; N_max = 20; d = size(x0,1);
	for k = 1:n
	    x(:,k+1) = x(:,k);
	    for j = 1:N_max
	        % Newton step 
	        F = x(:,k+1) - (x(:,k) + tau * feval(pdl,0,(x(:,k)+x(:,k+1))/2));
	        DF = eye(d) - tau * feval(dpdl,0,x(:,k+1));
	        dx = DF\F;
	        x(:,k+1) = x(:,k+1) - dx;
	        if (max(abs(dx)) < TOL) break; end;
	    end;
	end;
	h=plot(x(1,:),x(2,:),'k.-'); set(h,'linewidth',2,'markersize',20);