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

