Here it is.
Davide
function [tout,xout] = ode23s(FUN,tspan,x0,options)
% This file is intended for use with Octave.
% ode23s.m is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% ode23.m is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------
if nargin < 4, options=odeset();end
% Initialization
d=1/(2+ sqrt(2));
a=1/2;
e32=6+sqrt(2);
if size(options.RelTol)~=size([]) %user-defined relative tolerance
tol=options.RelTol;
else
tol = 1.e-3;
end
t = tspan(1);
tfinal = tspan(2);
if size(options.MaxStep)~=size([]) %user-defined max step size
hmax=options.MaxStep;
else
hmax = (tfinal - t)/2.5;
end
hmin = 16*eps*(tfinal - t);
if size(options.InitialStep)~=size([]) %user-defined initial step size
h=options.InitialStep;
else
h = (tfinal - t)/200; % initial guess at a step size
end
x = x0(:); % this always creates a column vector, x
tout = t; % first output time
xout = x.'; % first output solution
% The main loop
while (t < tfinal) && (h >= hmin)
if t + h > tfinal, h = tfinal - t; end
% Jacobian matrix, dfxpdp
J=dfxpdp(t,x,FUN);
T=(feval(FUN,x,t+hmin)-feval(FUN,x,t))/hmin;
% Wolfbrandt coefficient
W=eye(length(x0))-h*d*J;
% compute the slopes
F(:,1)=feval(FUN,x,t);
k(:,1)=W\(F(:,1)+ h*d*T);
F(:,2)=feval(FUN, x+a*h*k(:,1),t+a*h);
k(:,2)=W\((F(:,2) - k(:,1))) + k(:,1);
% compute the 2nd order estimate
x2=x + h*k(:,2);
% 3rd order, needed in error forumula
F(:,3)=feval(FUN,x2,t+h);
k(:,3)=W\(F(:,3)-e32*(k(:,2)-F(:,2))-2*(k(:,1)-F(:,1))+h*d*T);
% estimate the error
error = (h/6)*(norm(k(:,1)-2*k(:,2)+k(:,3)));
% Estimate the acceptable error
tau = tol*max(norm(x,'inf'),1.0);
% Update the solution only if the error is acceptable
if error <= tau
t = t + h;
x = x2; %no local extrapolation, FSAL
tout = [tout; t];
if size(options.Mass)~=size([]) %user-defined mass matrix
M=options.Mass;
xout = [xout;(M\x).'];
else
xout = [xout; x.'];
end
% Update the step size
if error == 0.0
error = 1e-16;
end
h = min(hmax, h*1.25); %auto-adaptive step update
else
h = max(hmin, h*0.5); %auto-adaptive step update
end
end;
if (t < tfinal)
disp('Step size grew too small.')
t, h, x
end------------------------------------------------------------------------------
Live Security Virtual Conference
Exclusive live event will cover all the ways today's security and
threat landscape has changed and how IT managers can respond. Discussions
will include endpoint security, mobile security and the latest in malware
threats. http://www.accelacomm.com/jaw/sfrnl04242012/114/50122263/
_______________________________________________
Octave-dev mailing list
Octave-dev@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/octave-dev
