function varargout = ode45(ode,tspan,y0,options,varargin)
solver_name = 'ode45';
if nargin < 4
options = [];
if nargin < 3
y0 = [];
if nargin < 2
tspan = [];
if nargin < 1
error('MATLAB
de45:NotEnoughInputs',...
'Not enough input arguments. See ODE45.');
end
end
end
end
% Stats
nsteps = 0;
nfailed = 0;
nfevals = 0;
% Output
FcnHandlesUsed = isa(ode,'function_handle');
output_sol = (FcnHandlesUsed && (nargout==1)); % sol = odeXX(...)
output_ty = (~output_sol && (nargout > 0)); % [t,y,...] = odeXX(...)
% There might be no output requested...
sol = []; f3d = [];
if output_sol
sol.solver = solver_name;
sol.extdata.odefun = ode;
sol.extdata.options = options;
sol.extdata.varargin = varargin;
end
% Handle solver arguments
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
options, threshold, rtol, normcontrol, normy, hmax, htry, htspan, dataType] = ...
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
nfevals = nfevals + 1;
% Handle the output
if nargout > 0
outputFcn = odeget(options,'OutputFcn',[],'fast');
else
outputFcn = odeget(options,'OutputFcn',@odeplot,'fast');
end
outputArgs = {};
if isempty(outputFcn)
haveOutputFcn = false;
else
haveOutputFcn = true;
outputs = odeget(options,'OutputSel',1:neq,'fast');
if isa(outputFcn,'function_handle')
% With MATLAB 6 syntax pass additional input arguments to outputFcn.
outputArgs = varargin;
end
end
refine = max(1,odeget(options,'Refine',4,'fast'));
if ntspan > 2
outputAt = 'RequestedPoints'; % output only at tspan points
elseif refine <= 1
outputAt = 'SolverSteps'; % computed points, no refinement
else
outputAt = 'RefinedSteps'; % computed points, with refinement
S = (1:refine-1) / refine;
end
printstats = strcmp(odeget(options,'Stats','off','fast'),'on');
% Handle the event function
[haveEventFcn,eventFcn,eventArgs,valt,teout,yeout,ieout] = ...
odeevents(FcnHandlesUsed,odeFcn,t0,y0,options,varargin);
% Handle the mass matrix
[Mtype, Mfun, Margs, M] = odemass(FcnHandlesUsed,odeFcn,t0,y0,options,varargin);
if Mtype>0 % non-trivial mass matrix
Msingular = odeget(options,'MassSingular','no','fast');
if strcmp(Msingular,'maybe')
warning('MATLAB
de45:MassSingularAssumedNo',['ODE45 assumes ' ...
'MassSingular is ''no''. See ODE15S or ODE23T.']);
elseif strcmp(Msingular,'yes')
error('MATLAB
de45:MassSingularYes',...
['MassSingular cannot be ''yes'' for this solver. See ODE15S '...
' or ODE23T.']);
end
if Mtype == 1
[L,U] = lu(M);
else
L = [];
U = [];
end
% Incorporate the mass matrix into odeFcn and odeArgs.
[odeFcn,odeArgs] = odemassexplicit(FcnHandlesUsed,Mtype,odeFcn,odeArgs,Mfun,Margs,L,U);
f0 = feval(odeFcn,t0,y0,odeArgs{:});
nfevals = nfevals + 1;
end
% Non-negative solution components
idxNonNegative = odeget(options,'NonNegative',[],'fast');
nonNegative = ~isempty(idxNonNegative);
if nonNegative % modify the derivative function
[odeFcn,thresholdNonNegative] = odenonnegative(odeFcn,y0,threshold,idxNonNegative);
f0 = feval(odeFcn,t0,y0,odeArgs{:});
nfevals = nfevals + 1;
end
t = t0;
y = y0;
% Allocate memory if we're generating output.
nout = 0;
tout = []; yout = [];
if nargout > 0
if output_sol
chunk = min(max(100,50*refine), refine+floor((2^11)/neq));
tout = zeros(1,chunk,dataType);
yout = zeros(neq,chunk,dataType);
f3d = zeros(neq,7,chunk,dataType);
else
if ntspan > 2 % output only at tspan points
tout = zeros(1,ntspan,dataType);
yout = zeros(neq,ntspan,dataType);
else % alloc in chunks
chunk = min(max(100,50*refine), refine+floor((2^13)/neq));
tout = zeros(1,chunk,dataType);
yout = zeros(neq,chunk,dataType);
end
end
nout = 1;
tout(nout) = t;
yout(:,nout) = y;
end
% Initialize method parameters.
pow = 1/5;
A = [1/5, 3/10, 4/5, 8/9, 1, 1];
B = [
1/5 3/40 44/45 19372/6561 9017/3168 35/384
0 9/40 -56/15 -25360/2187 -355/33 0
0 0 32/9 64448/6561 46732/5247 500/1113
0 0 0 -212/729 49/176 125/192
0 0 0 0 -5103/18656 -2187/6784
0 0 0 0 0 11/84
0 0 0 0 0 0
];
E = [71/57600; 0; -71/16695; 71/1920; -17253/339200; 22/525; -1/40];
f = zeros(neq,7,dataType);
hmin = 16*eps(t);
if isempty(htry)
% Compute an initial step size h using y'(t).
absh = min(hmax, htspan);
if normcontrol
rh = (norm(f0) / max(normy,threshold)) / (0.8 * rtol^pow);
else
rh = norm(f0 ./ max(abs(y),threshold),inf) / (0.8 * rtol^pow);
end
if absh * rh > 1
absh = 1 / rh;
end
absh = max(absh, hmin);
else
absh = min(hmax, max(hmin, htry));
end
f(:,1) = f0;
% Initialize the output function.
if haveOutputFcn
feval(outputFcn,[t tfinal],y(outputs),'init',outputArgs{:});
end
% THE MAIN LOOP
done = false;
while ~done
% By default, hmin is a small number such that t+hmin is only slightly
% different than t. It might be 0 if t is 0.
hmin = 16*eps(t);
absh = min(hmax, max(hmin, absh)); % couldn't limit absh until new hmin
h = tdir * absh;
% Stretch the step if within 10% of tfinal-t.
if 1.1*absh >= abs(tfinal - t)
h = tfinal - t;
absh = abs(h);
done = true;
end
% LOOP FOR ADVANCING ONE STEP.
nofailed = true; % no failed attempts
while true
hA = h * A;
hB = h * B;
f(:,2) = feval(odeFcn,t+hA(1),y+f*hB(:,1),odeArgs{:});
f(:,3) = feval(odeFcn,t+hA(2),y+f*hB(:,2),odeArgs{:});
f(:,4) = feval(odeFcn,t+hA(3),y+f*hB(:,3),odeArgs{:});
f(:,5) = feval(odeFcn,t+hA(4),y+f*hB(:,4),odeArgs{:});
f(:,6) = feval(odeFcn,t+hA(5),y+f*hB(:,5),odeArgs{:});
tnew = t + hA(6);
if done
tnew = tfinal; % Hit end point exactly.
end
h = tnew - t; % Purify h.
ynew = y + f*hB(:,6);
f(:,7) = feval(odeFcn,tnew,ynew,odeArgs{:});
nfevals = nfevals + 6;
% Estimate the error.
NNrejectStep = false;
if normcontrol
normynew = norm(ynew);
errwt = max(max(normy,normynew),threshold);
err = absh * (norm(f * E) / errwt);
if nonNegative && (err <= rtol) && any(ynew(idxNonNegative)<0)
errNN = norm( max(0,-ynew(idxNonNegative)) ) / errwt ;
if errNN > rtol
err = errNN;
NNrejectStep = true;
end
end
else
err = absh * norm((f * E) ./ max(max(abs(y),abs(ynew)),threshold),inf);
if nonNegative && (err <= rtol) && any(ynew(idxNonNegative)<0)
errNN = norm( max(0,-ynew(idxNonNegative)) ./ thresholdNonNegative, inf);
if errNN > rtol
err = errNN;
NNrejectStep = true;
end
end
end
if err > rtol % Failed step
nfailed = nfailed + 1;
if absh <= hmin
warning('MATLAB
de45:IntegrationTolNotMet',['Failure at t=%e. ' ...
'Unable to meet integration tolerances without reducing ' ...
'the step size below the smallest value allowed (%e) ' ...
'at time t.'],t,hmin);
solver_output = odefinalize(solver_name, sol,...
outputFcn, outputArgs,...
printstats, [nsteps, nfailed, nfevals],...
nout, tout, yout,...
haveEventFcn, teout, yeout, ieout,...
{f3d,idxNonNegative});
if nargout > 0
varargout = solver_output;
end
return;
end
if nofailed
nofailed = false;
if NNrejectStep
absh = max(hmin, 0.5*absh);
else
absh = max(hmin, absh * max(0.1, 0.8*(rtol/err)^pow));
end
else
absh = max(hmin, 0.5 * absh);
end
h = tdir * absh;
done = false;
else % Successful step
NNreset_f7 = false;
if nonNegative && any(ynew(idxNonNegative)<0)
ynew(idxNonNegative) = max(ynew(idxNonNegative),0);
if normcontrol
normynew = norm(ynew);
end
NNreset_f7 = true;
end
break;
end
end
nsteps = nsteps + 1;
if haveEventFcn
[te,ye,ie,valt,stop] = ...
odezero(@ntrp45,eventFcn,eventArgs,valt,t,y,tnew,ynew,t0,h,f,idxNonNegative);
if ~isempty(te)
if output_sol || (nargout > 2)
teout = [teout, te];
yeout = [yeout, ye];
ieout = [ieout, ie];
end
if stop % Stop on a terminal event.
% Adjust the interpolation data to [t te(end)].
% Update the derivatives using the interpolating polynomial.
taux = t + (te(end) - t)*A;
[ignore,f(:,2:7)] = ntrp45(taux,t,y,[],[],h,f,idxNonNegative);
tnew = te(end);
ynew = ye(:,end);
h = tnew - t;
done = true;
end
end
end
if output_sol
nout = nout + 1;
if nout > length(tout)
tout = [tout, zeros(1,chunk,dataType)]; % requires chunk >= refine
yout = [yout, zeros(neq,chunk,dataType)];
f3d = cat(3,f3d,zeros(neq,7,chunk,dataType));
end
tout(nout) = tnew;
yout(:,nout) = ynew;
f3d(:,:,nout) = f;
end
if output_ty || haveOutputFcn
switch outputAt
case 'SolverSteps' % computed points, no refinement
nout_new = 1;
tout_new = tnew;
yout_new = ynew;
case 'RefinedSteps' % computed points, with refinement
tref = t + (tnew-t)*S;
nout_new = refine;
tout_new = [tref, tnew];
yout_new = [ntrp45(tref,t,y,[],[],h,f,idxNonNegative), ynew];
case 'RequestedPoints' % output only at tspan points
nout_new = 0;
tout_new = [];
yout_new = [];
while next <= ntspan
if tdir * (tnew - tspan(next)) < 0
if haveEventFcn && stop % output tstop,ystop
nout_new = nout_new + 1;
tout_new = [tout_new, tnew];
yout_new = [yout_new, ynew];
end
break;
end
nout_new = nout_new + 1;
tout_new = [tout_new, tspan(next)];
if tspan(next) == tnew
yout_new = [yout_new, ynew];
else
yout_new = [yout_new, ntrp45(tspan(next),t,y,[],[],h,f,idxNonNegative)];
end
next = next + 1;
end
end
if nout_new > 0
if output_ty
oldnout = nout;
nout = nout + nout_new;
if nout > length(tout)
tout = [tout, zeros(1,chunk,dataType)]; % requires chunk >= refine
yout = [yout, zeros(neq,chunk,dataType)];
end
idx = oldnout+1:nout;
tout(idx) = tout_new;
yout(:,idx) = yout_new;
end
if haveOutputFcn
stop = feval(outputFcn,tout_new,yout_new(outputs,
,'',outputArgs{:});
if stop
done = true;
end
end
end
end
if done
break
end
% If there were no failures compute a new h.
if nofailed
% Note that absh may shrink by 0.8, and that err may be 0.
temp = 1.25*(err/rtol)^pow;
if temp > 0.2
absh = absh / temp;
else
absh = 5.0*absh;
end
end
% Advance the integration one step.
t = tnew;
y = ynew;
if normcontrol
normy = normynew;
end
if NNreset_f7
% Used f7 for unperturbed solution to interpolate.
% Now reset f7 to move along constraint.
f(:,7) = feval(odeFcn,tnew,ynew,odeArgs{:});
nfevals = nfevals + 1;
end
f(:,1) = f(:,7); % Already have f(tnew,ynew)
end
solver_output = odefinalize(solver_name, sol,...
outputFcn, outputArgs,...
printstats, [nsteps, nfailed, nfevals],...
nout, tout, yout,...
haveEventFcn, teout, yeout, ieout,...
{f3d,idxNonNegative});
if nargout > 0
varargout = solver_output;
end