function [tout, yout] = ppsolvep(FunFcn, y0, window,magn,sign,steps,tol,itlim)

% PPSOLVEP	Solve differential equations, higher order method.
%	ODE45 integrates a system of ordinary differential equations using
%	4th and 5th order Runge-Kutta formulas.
%
%	See also ODE23, ODEDEMO.

%	C.B. Moler, 3-25-87.
%	Copyright (c) 1984-92 by The MathWorks, Inc.

% The Fehlberg coefficients:

% alpha = [1/4  3/8  12/13  1  1/2]';   Not needed for autonomous systems.

beta  = [ [    1      0      0     0      0    0]/4
          [    3      9      0     0      0    0]/32
          [ 1932  -7200   7296     0      0    0]/2197
          [ 8341 -32832  29440  -845      0    0]/4104
          [-6080  41040 -28352  9295  -5643    0]/20520 ]';
gamma = [ [902880  0  3953664  3855735  -1371249  277020]/7618050
          [ -2090  0    22528    21970    -15048  -27360]/752400 ]';
pow = 1/5;

if nargin < 7, tol = 1e-6; end

% itlim=1000;

hmin = 1e-10;
% NN = steps;		% The minimum number of trajectory elements

window = window(:);
dwindow = [window(1), window(3), -window(2), -window(4)].';

DY=[window(2)-window(1);window(4)-window(3)]; % = [dx,dy]

cwindow = dwindow - magn*[DY;DY];
blocksize = 100;

index = 0;
lastindex = blocksize + 1;

yy0=[y0(:);-y0(:)];

inside = yy0-dwindow;ins = all(inside>=0);

% Initialization
t = 0;
h = 0.00001*sign;
y = y0(:);
f = y*zeros(1,6);
tout = [t;zeros(blocksize,1)];
yout = [y.';zeros(blocksize,2)];

tau = tol * max(norm(y, 'inf'), 1);

% Set up the plot routine.

X=y(1);Y=y(2);
ppp = plot(X,Y,'-r','Erasemode','none','clip','on');

% Initialize the flags.

windowflag = all(yy0 -cwindow >=0);
sinkflag = 1;
orbitflag = 1;
stepflag = 1;

theta = rand*2*pi;
R=[cos(theta),sin(theta);-sin(theta),cos(theta)];
turn=[];
turnk=0;
vv=R*DY;

% The closeness required perpendicular to a closed orbit.
eps1=1/(norm(R(:,1)./DY)*1000);
	
% The closeness required parallel to the orbit.
eps2=2/(norm(R(:,2)./DY)*10);

% The closeness required to detect an equilibrium point.
eps3=1/(300);

N=0;			

% The main loop
while ((windowflag)&(sinkflag)&(orbitflag)&(stepflag)&(itlim))

	% Compute the slopes
	temp = feval(FunFcn,y);
	f(:,1) = temp(:);
	for j = 1:5
	 	temp = feval(FunFcn, y+h*f*beta(:,j));
	 	f(:,j+1) = temp(:);
	end
	
	% Estimate the error and the acceptable error
	delta = norm(h*f*gamma(:,2));
	tau = tol*max(norm(y),1.0);
	itlim = itlim-1;
	Mf=f*gamma(:,1);
	MMf=norm(Mf./DY);
	
% Update the solution only if the error is acceptable.

	if ( (delta <= tau) )
	 	tn = t + h;
	 	yn = y + h*Mf;
		index = index +1;
		if ( index > lastindex)
			tout = [tout;zeros(blocksize,1)];
			yout = [yout;zeros(blocksize,2)];
			lastindex = lastindex + blocksize;
		end

        tout(index) = tn;
        yout(index,:) = yn.';

	 % Update the plot if either y or yn is in the display window.
	 % If only one is in the window, it will be necessary to clip
	 % the plot.
	 
	 	insiden = [yn(:);-yn(:)] - dwindow;insn = all(insiden>=0);
		if (ins & insn)
			set(ppp,'Xdata',[y(1),yn(1)],'Ydata',[y(2),yn(2)]);
		 	drawnow
		elseif (ins)
			kkk = find(insiden<0);
			lll = -inside(kkk)./(insiden(kkk)-inside(kkk));
			lll = min(lll);
			ynb = y + lll*(yn-y);
			set(ppp,'Xdata',[y(1),ynb(1)],'Ydata',[y(2),ynb(2)]);
			drawnow
		elseif (insn)
			kkk = find(inside<0);
			lll = -insiden(kkk)./(inside(kkk)-insiden(kkk));
			lll = min(lll);
			ynb = yn + lll*(y-yn);
			set(ppp,'Xdata',[ynb(1),yn(1)],'Ydata',[ynb(2),yn(2)]);
		 	drawnow
		end
			
% Update the flags
    
		windowflag = all([yn(:);-yn(:)] - cwindow >=0);

      	if (index > 30)
			sinkflag= (MMf>eps3);	% (norm((y-yn)./DY)>eps3);
	 	end
		
		if (index == 1)
			pp=R*yn;
		elseif (index == 2)
			qq=R*yn;
		else
			rr = R*yn;
			if( (pp(1)<qq(1)) & (rr(1)<qq(1)) )
				kk=0;
				while ( (kk<turnk) & (orbitflag) )
					kk = kk+1;
					if ( (abs(qq(1)-turn(1,kk))<eps1) &...
						(abs(qq(2)-turn(2,kk))<eps2) )
						
						orbitflag = 0;
					end
				end
				if (orbitflag)
					turn = [turn,qq];
					turnk = turnk + 1;
				end
			end
			pp = qq;
			qq = rr;
		end
	
	 	t=tn;y=yn;
	 	inside = insiden;ins = insn;
	end

      % Update the step size

    if (delta ~= 0.0)
        absh = min(0.8*abs(h)*(tau/delta)^pow, 1);  % 1/(NN*MMf+eps));
	 	h=sign*absh;
	 	stepflag = (absh>=hmin);
    else
	 	h= sign;   % /(NN*MMf+eps);
	end
end;


%if (windowflag == 0)
%	disp('The orbit has left the computation window.');
%end

if (sinkflag == 0)
	disp(['The orbit ends in a possible equilibrium point near (',...
			num2str(yn(1)), ', ', num2str(yn(2)),').']);
end

if (orbitflag == 0)
	disp('A nearly closed orbit was detected.');
end

if (itlim<=0)
	disp(['Maximum number of iterations reached.']);
end

if (stepflag == 0)
	disp('A step size smaller than the minimum required.');
end

tout = tout(1:index);
yout = yout(1:index,:);