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

% usage: [tout, yout] = dfsolve(FunFcn, y0, window,magn,sign,steps,tol,itlim)
%
% DFSOLVE	Solve differential equations, higher order method.
%	DFSOLVE integrates a single ordinary differential equation 
%	using 4th and 5th order Runge-Kutta formulas.  This routine 
%	is designed to be used only with DFIELD.
%

%	Based on ODE$%.m written by
%	C.B. Moler, 3-25-87.
%	Copyright (c) 1984-92 by The MathWorks, Inc.

%	Last modified by John C. Polking, Rice Unviersity
%	February 3, 1995

% The Fehlberg coefficients:

alpha = [1/4  3/8  12/13  1  1/2]';

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;

% 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;

hmin = DY(1)*1e-7;
hmax = DY(1)/20;

index = 0;
lastindex = blocksize + 1;

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

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

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

tau = tol * max(abs(y), 1);

% Set up the plot routine.


ppp = plot(t,y,'-y','Erasemode','none','clip','on');

% Initialize the flags.

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

N=0;			

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

	% Compute the slopes
	temp = feval(FunFcn,t,y);
	f(1) = temp;
	for j = 1:5
	 	temp = feval(FunFcn, t+alpha(j)*h, y+h*f*beta(:,j));
	 	f(j+1) = temp;
	end
	
	% Estimate the error and the acceptable error
	delta = abs(h*f*gamma(:,2));
	tau = tol*max(abs(y),1.0);
	itlim = itlim-1;
	Mf=f*gamma(:,1);
	MMf=abs(Mf/DY(2));
	
% 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,1)];
			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 = [tn;yn;-tn;-yn] - dwindow;insn = all(insiden>=0);
		if (ins & insn)
			set(ppp,'Xdata',[t,tn],'Ydata',[y,yn]);
		 	drawnow
		elseif (ins)
			kkk = find(insiden<0);
			lll = -inside(kkk)./(insiden(kkk)-inside(kkk));
			lll = min(lll);
			ynb = y + lll*(yn-y);
			tnb = t + lll*(tn-t);
			set(ppp,'Xdata',[t,tnb],'Ydata',[y,ynb]);
			drawnow
		elseif (insn)
			kkk = find(inside<0);
			lll = -insiden(kkk)./(inside(kkk)-insiden(kkk));
			lll = min(lll);
			ynb = yn + lll*(y-yn);
			tnb = tn + lll*(t-tn);
			set(ppp,'Xdata',[tnb,tn],'Ydata',[ynb,yn]);
		 	drawnow
		end
			
% Update the window flag
    
		windowflag = all([tn;yn;-tn;-yn] - cwindow >=0);
	
	 	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, hmax);  % 1/(NN*MMf+eps));
	 	h=sign*absh;
	 	stepflag = (absh>=hmin);
    else
	 	h= sign*hmax;  
	end
end;


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

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

if (stepflag == 0)
	disp('A step size smaller than the minimum required.');
	disp(['Singularity possible near (',...
		num2str(tn), ', ', num2str(yn),').'])
end

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