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