# HG changeset patch
# User Doug Stewart
-%# OdePkg - A package for solving ordinary differential equations and more
-%#
-%# This program 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 of the License, or
-%# (at your option) any later version.
-%#
-%# This program 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.
-%#
-%# You should have received a copy of the GNU General Public License
-%# along with this program; If not, see .
+## Copyright (C) 2006-2012, Thomas Treichl
+## OdePkg - A package for solving ordinary differential equations and more
+##
+## This program 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 of the License, or
+## (at your option) any later version.
+##
+## This program 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.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; If not, see .
-%# -*- texinfo -*-
-%# @deftypefn {Function File} address@hidden =} ode23 (@var{@@fun}, @var{slot}, @var{init}, address@hidden, address@hidden, @var{par2}, @dots{}])
-%# @deftypefnx {Command} address@hidden =} ode23 (@var{@@fun}, @var{slot}, @var{init}, address@hidden, address@hidden, @var{par2}, @dots{}])
-%# @deftypefnx {Command} address@hidden, @var{y}, address@hidden, @var{ye}, @var{ie}]] =} ode23 (@var{@@fun}, @var{slot}, @var{init}, address@hidden, address@hidden, @var{par2}, @dots{}])
-%#
-%# This function file can be used to solve a set of non--stiff ordinary differential equations (non--stiff ODEs) or non--stiff differential algebraic equations (non--stiff DAEs) with the well known explicit Runge--Kutta method of order (2,3).
-%#
-%# If this function is called with no return argument then plot the solution over time in a figure window while solving the set of ODEs that are defined in a function and specified by the function handle @var{@@fun}. The second input argument @var{slot} is a double vector that defines the time slot, @var{init} is a double vector that defines the initial values of the states, @var{opt} can optionally be a structure array that keeps the options created with the command @command{odeset} and @var{par1}, @var{par2}, @dots{} can optionally be other input arguments of any type that have to be passed to the function defined by @var{@@fun}.
-%#
-%# If this function is called with one return argument then return the solution @var{sol} of type structure array after solving the set of ODEs. The solution @var{sol} has the fields @var{x} of type double column vector for the steps chosen by the solver, @var{y} of type double column vector for the solutions at each time step of @var{x}, @var{solver} of type string for the solver name and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector that keep the informations of the event function if an event function handle is set in the option argument @var{opt}.
-%#
-%# If this function is called with more than one return argument then return the time stamps @var{t}, the solution values @var{y} and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector.
-%#
-%# For example, solve an anonymous implementation of the Van der Pol equation
-%#
-%# @example
-%# fvdb = @@(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
-%#
-%# vopt = odeset ("RelTol", 1e-3, "AbsTol", 1e-3, \
-%# "NormControl", "on", "OutputFcn", @@odeplot);
-%# ode23 (fvdb, [0 20], [2 0], vopt);
-%# @end example
-%# @end deftypefn
-%#
-%# @seealso{odepkg}
+## -*- texinfo -*-
+## @deftypefn {Function File} address@hidden =} ode23 (@var{@@fun}, @var{slot}, @var{init}, address@hidden, address@hidden, @var{par2}, @dots{}])
+## @deftypefnx {Command} address@hidden =} ode23 (@var{@@fun}, @var{slot}, @var{init}, address@hidden, address@hidden, @var{par2}, @dots{}])
+## @deftypefnx {Command} address@hidden, @var{y}, address@hidden, @var{ye}, @var{ie}]] =} ode23 (@var{@@fun}, @var{slot}, @var{init}, address@hidden, address@hidden, @var{par2}, @dots{}])
+##
+## This function file can be used to solve a set of non--stiff ordinary differential equations (non--stiff ODEs) or non--stiff differential algebraic equations (non--stiff DAEs) with the well known explicit Runge--Kutta method of order (2,3).
+##
+## If this function is called with no return argument then plot the solution over time in a figure window while solving the set of ODEs that are defined in a function and specified by the function handle @var{@@fun}. The second input argument @var{slot} is a double vector that defines the time slot, @var{init} is a double vector that defines the initial values of the states, @var{opt} can optionally be a structure array that keeps the options created with the command @command{odeset} and @var{par1}, @var{par2}, @dots{} can optionally be other input arguments of any type that have to be passed to the function defined by @var{@@fun}.
+##
+## If this function is called with one return argument then return the solution @var{sol} of type structure array after solving the set of ODEs. The solution @var{sol} has the fields @var{x} of type double column vector for the steps chosen by the solver, @var{y} of type double column vector for the solutions at each time step of @var{x}, @var{solver} of type string for the solver name and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector that keep the informations of the event function if an event function handle is set in the option argument @var{opt}.
+##
+## If this function is called with more than one return argument then return the time stamps @var{t}, the solution values @var{y} and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector.
+##
+## For example, solve an anonymous implementation of the Van der Pol equation
+##
+## @example
+## fvdb = @@(vt, vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
+##
+## vopt = odeset ("RelTol", 1e-3, "AbsTol", 1e-3, ...
+## "NormControl", "on", "OutputFcn", @@odeplot);
+## ode23 (fvdb, [0 20], [2 0], vopt);
+## @end example
+## @end deftypefn
+##
+## @seealso{odepkg}
-%# ChangeLog:
-%# 20010703 the function file "ode23.m" was written by Marc Compere
-%# under the GPL for the use with this software. This function has been
-%# taken as a base for the following implementation.
-%# 20060810, Thomas Treichl
-%# This function was adapted to the new syntax that is used by the
-%# new OdePkg for Octave and is compatible to Matlab's ode23.
+## ChangeLog:
+## 20010703 the function file "ode23.m" was written by Marc Compere
+## under the GPL for the use with this software. This function has been
+## taken as a base for the following implementation.
+## 20060810, Thomas Treichl
+## This function was adapted to the new syntax that is used by the
+## new OdePkg for Octave and is compatible to Matlab's ode23.
function [varargout] = ode23 (vfun, vslot, vinit, varargin)
- if (nargin == 0) %# Check number and types of all input arguments
- help ('ode23');
- error ('OdePkg:InvalidArgument', ...
- 'Number of input arguments must be greater than zero');
+ if (nargin == 0) # Check number and types of all input arguments
+ help ("ode23");
+ error ("OdePkg:InvalidArgument", ...
+ "Number of input arguments must be greater than zero");
elseif (nargin < 3)
print_usage;
- elseif ~(isa (vfun, 'function_handle') || isa (vfun, 'inline'))
- error ('OdePkg:InvalidArgument', ...
- 'First input argument must be a valid function handle');
+ elseif ! (isa (vfun, "function_handle") || isa (vfun, "inline"))
+ error ("OdePkg:InvalidArgument", ...
+ "First input argument must be a valid function handle");
- elseif (~isvector (vslot) || length (vslot) < 2)
- error ('OdePkg:InvalidArgument', ...
- 'Second input argument must be a valid vector');
+ elseif (! isvector (vslot) || length (vslot) < 2)
+ error ("OdePkg:InvalidArgument", ...
+ "Second input argument must be a valid vector");
- elseif (~isvector (vinit) || ~isnumeric (vinit))
- error ('OdePkg:InvalidArgument', ...
- 'Third input argument must be a valid numerical value');
+ elseif (! isvector (vinit) || ! isnumeric (vinit))
+ error ("OdePkg:InvalidArgument", ...
+ "Third input argument must be a valid numerical value");
elseif (nargin >= 4)
- if (~isstruct (varargin{1}))
- %# varargin{1:len} are parameters for vfun
+ if (! isstruct (varargin{1}))
+ # varargin{1:len} are parameters for vfun
vodeoptions = odeset;
vfunarguments = varargin;
elseif (length (varargin) > 1)
- %# varargin{1} is an OdePkg options structure vopt
- vodeoptions = odepkg_structure_check (varargin{1}, 'ode23');
+ # varargin{1} is an OdePkg options structure vopt
+ vodeoptions = odepkg_structure_check (varargin{1}, "ode23");
vfunarguments = {varargin{2:length(varargin)}};
- else %# if (isstruct (varargin{1}))
- vodeoptions = odepkg_structure_check (varargin{1}, 'ode23');
+ else # if (isstruct (varargin{1}))
+ vodeoptions = odepkg_structure_check (varargin{1}, "ode23");
vfunarguments = {};
- end
+ endif
- else %# if (nargin == 3)
+ else # if (nargin == 3)
vodeoptions = odeset;
vfunarguments = {};
- end
+ endif
- %# Start preprocessing, have a look which options are set in
- %# vodeoptions, check if an invalid or unused option is set
- vslot = vslot(:).'; %# Create a row vector
- vinit = vinit(:).'; %# Create a row vector
- if (length (vslot) > 2) %# Step size checking
+ ## Start preprocessing, have a look which options are set in
+ ## vodeoptions, check if an invalid or unused option is set
+ vslot = vslot(:).'; # Create a row vector
+ vinit = vinit(:).'; # Create a row vector
+ if (length (vslot) > 2) # Step size checking
vstepsizefixed = true;
else
vstepsizefixed = false;
- end
+ endif
- %# Get the default options that can be set with 'odeset' temporarily
+ ## Get the default options that can be set with "odeset" temporarily
vodetemp = odeset;
- %# Implementation of the option RelTol has been finished. This option
- %# can be set by the user to another value than default value.
- if (isempty (vodeoptions.RelTol) && ~vstepsizefixed)
+ ## Implementation of the option RelTol has been finished. This option
+ ## can be set by the user to another value than default value.
+ if (isempty (vodeoptions.RelTol) && ! vstepsizefixed)
vodeoptions.RelTol = 1e-6;
- warning ('OdePkg:InvalidArgument', ...
- 'Option "RelTol" not set, new value %f is used', vodeoptions.RelTol);
- elseif (~isempty (vodeoptions.RelTol) && vstepsizefixed)
- warning ('OdePkg:InvalidArgument', ...
- 'Option "RelTol" will be ignored if fixed time stamps are given');
- end
+ warning ("OdePkg:InvalidArgument", ...
+ "Option 'RelTol' not set, new value %f is used", vodeoptions.RelTol);
+ elseif (! isempty (vodeoptions.RelTol) && vstepsizefixed)
+ warning ("OdePkg:InvalidArgument", ...
+ "Option 'RelTol' will be ignored if fixed time stamps are given");
+ endif
- %# Implementation of the option AbsTol has been finished. This option
- %# can be set by the user to another value than default value.
- if (isempty (vodeoptions.AbsTol) && ~vstepsizefixed)
+ ## Implementation of the option AbsTol has been finished. This option
+ # can be set by the user to another value than default value.
+ if (isempty (vodeoptions.AbsTol) && ! vstepsizefixed)
vodeoptions.AbsTol = 1e-6;
- warning ('OdePkg:InvalidArgument', ...
- 'Option "AbsTol" not set, new value %f is used', vodeoptions.AbsTol);
- elseif (~isempty (vodeoptions.AbsTol) && vstepsizefixed)
- warning ('OdePkg:InvalidArgument', ...
- 'Option "AbsTol" will be ignored if fixed time stamps are given');
+ warning ("OdePkg:InvalidArgument", ...
+ "Option 'AbsTol' not set, new value %f is used", vodeoptions.AbsTol);
+ elseif (! isempty (vodeoptions.AbsTol) && vstepsizefixed)
+ warning ("OdePkg:InvalidArgument", ...
+ "Option 'AbsTol' will be ignored if fixed time stamps are given");
else
- vodeoptions.AbsTol = vodeoptions.AbsTol(:); %# Create column vector
- end
+ vodeoptions.AbsTol = vodeoptions.AbsTol(:); # Create column vector
+ endif
- %# Implementation of the option NormControl has been finished. This
- %# option can be set by the user to another value than default value.
- if (strcmp (vodeoptions.NormControl, 'on')) vnormcontrol = true;
- else vnormcontrol = false; end
+ ## Implementation of the option NormControl has been finished. This
+ ## option can be set by the user to another value than default value.
+ if (strcmp (vodeoptions.NormControl, "on")) vnormcontrol = true;
+ else vnormcontrol = false;
+ endif
- %# Implementation of the option NonNegative has been finished. This
- %# option can be set by the user to another value than default value.
- if (~isempty (vodeoptions.NonNegative))
+ ## Implementation of the option NonNegative has been finished. This
+ ## option can be set by the user to another value than default value.
+ if (! isempty (vodeoptions.NonNegative))
if (isempty (vodeoptions.Mass)), vhavenonnegative = true;
else
vhavenonnegative = false;
- warning ('OdePkg:InvalidArgument', ...
- 'Option "NonNegative" will be ignored if mass matrix is set');
- end
+ warning ("OdePkg:InvalidArgument", ...
+ "Option 'NonNegative' will be ignored if mass matrix is set");
+ endif
else vhavenonnegative = false;
- end
+ endif
- %# Implementation of the option OutputFcn has been finished. This
- %# option can be set by the user to another value than default value.
+ ## Implementation of the option OutputFcn has been finished. This
+ ## option can be set by the user to another value than default value.
if (isempty (vodeoptions.OutputFcn) && nargout == 0)
vodeoptions.OutputFcn = @odeplot;
vhaveoutputfunction = true;
elseif (isempty (vodeoptions.OutputFcn)), vhaveoutputfunction = false;
else vhaveoutputfunction = true;
- end
+ endif
- %# Implementation of the option OutputSel has been finished. This
- %# option can be set by the user to another value than default value.
- if (~isempty (vodeoptions.OutputSel)), vhaveoutputselection = true;
- else vhaveoutputselection = false; end
+ ## Implementation of the option OutputSel has been finished. This
+ ## option can be set by the user to another value than default value.
+ if (! isempty (vodeoptions.OutputSel)), vhaveoutputselection = true;
+ else vhaveoutputselection = false;
+ endif
- %# Implementation of the option OutputSave has been finished. This
- %# option can be set by the user to another value than default value.
+ ## Implementation of the option OutputSave has been finished. This
+ ## option can be set by the user to another value than default value.
if (isempty (vodeoptions.OutputSave)), vodeoptions.OutputSave = 1;
- end
+ endif
- %# Implementation of the option Refine has been finished. This option
- %# can be set by the user to another value than default value.
+ ## Implementation of the option Refine has been finished. This option
+ ## can be set by the user to another value than default value.
if (vodeoptions.Refine > 0), vhaverefine = true;
- else vhaverefine = false; end
+ else vhaverefine = false;
+ endif
- %# Implementation of the option Stats has been finished. This option
- %# can be set by the user to another value than default value.
+ ## Implementation of the option Stats has been finished. This option
+ ## can be set by the user to another value than default value.
- %# Implementation of the option InitialStep has been finished. This
- %# option can be set by the user to another value than default value.
- if (isempty (vodeoptions.InitialStep) && ~vstepsizefixed)
- vodeoptions.InitialStep = (vslot(1,2) - vslot(1,1)) / 10;
+ ## Implementation of the option InitialStep has been finished. This
+ ## option can be set by the user to another value than default value.
+ if (isempty (vodeoptions.InitialStep) && ! vstepsizefixed)
+ vodeoptions.InitialStep = (vslot(1, 2) - vslot(1, 1)) / 10;
vodeoptions.InitialStep = vodeoptions.InitialStep / 10^vodeoptions.Refine;
- warning ('OdePkg:InvalidArgument', ...
- 'Option "InitialStep" not set, new value %f is used', vodeoptions.InitialStep);
- end
+ warning ("OdePkg:InvalidArgument", ...
+ "Option 'InitialStep' not set, new value %f is used", vodeoptions.InitialStep);
+ endif
- %# Implementation of the option MaxStep has been finished. This option
- %# can be set by the user to another value than default value.
- if (isempty (vodeoptions.MaxStep) && ~vstepsizefixed)
- vodeoptions.MaxStep = abs (vslot(1,2) - vslot(1,1)) / 10;
- warning ('OdePkg:InvalidArgument', ...
- 'Option "MaxStep" not set, new value %f is used', vodeoptions.MaxStep);
- end
+ ## Implementation of the option MaxStep has been finished. This option
+ ## can be set by the user to another value than default value.
+ if (isempty (vodeoptions.MaxStep) && ! vstepsizefixed)
+ vodeoptions.MaxStep = abs (vslot(1, 2) - vslot(1, 1)) / 10;
+ warning ("OdePkg:InvalidArgument", ...
+ "Option 'MaxStep' not set, new value %f is used", vodeoptions.MaxStep);
+ endif
- %# Implementation of the option Events has been finished. This option
- %# can be set by the user to another value than default value.
- if (~isempty (vodeoptions.Events)), vhaveeventfunction = true;
- else vhaveeventfunction = false; end
+ ## Implementation of the option Events has been finished. This option
+ ## can be set by the user to another value than default value.
+ if (! isempty (vodeoptions.Events)), vhaveeventfunction = true;
+ else vhaveeventfunction = false;
+ endif
- %# The options 'Jacobian', 'JPattern' and 'Vectorized' will be ignored
- %# by this solver because this solver uses an explicit Runge-Kutta
- %# method and therefore no Jacobian calculation is necessary
- if (~isequal (vodeoptions.Jacobian, vodetemp.Jacobian))
- warning ('OdePkg:InvalidArgument', ...
- 'Option "Jacobian" will be ignored by this solver');
- end
- if (~isequal (vodeoptions.JPattern, vodetemp.JPattern))
- warning ('OdePkg:InvalidArgument', ...
- 'Option "JPattern" will be ignored by this solver');
- end
- if (~isequal (vodeoptions.Vectorized, vodetemp.Vectorized))
- warning ('OdePkg:InvalidArgument', ...
- 'Option "Vectorized" will be ignored by this solver');
- end
- if (~isequal (vodeoptions.NewtonTol, vodetemp.NewtonTol))
- warning ('OdePkg:InvalidArgument', ...
- 'Option "NewtonTol" will be ignored by this solver');
- end
- if (~isequal (vodeoptions.MaxNewtonIterations,...
+ ## The options "Jacobian", "JPattern" and "Vectorized" will be ignored
+ ## by this solver because this solver uses an explicit Runge-Kutta
+ ## method and therefore no Jacobian calculation is necessary
+ if (! isequal (vodeoptions.Jacobian, vodetemp.Jacobian))
+ warning ("OdePkg:InvalidArgument", ...
+ "Option 'Jacobian' will be ignored by this solver");
+ endif
+ if (! isequal (vodeoptions.JPattern, vodetemp.JPattern))
+ warning ("OdePkg:InvalidArgument", ...
+ "Option 'JPattern' will be ignored by this solver");
+ endif
+ if (! isequal (vodeoptions.Vectorized, vodetemp.Vectorized))
+ warning ("OdePkg:InvalidArgument", ...
+ "Option 'Vectorized' will be ignored by this solver");
+ endif
+ if (! isequal (vodeoptions.NewtonTol, vodetemp.NewtonTol))
+ warning ("OdePkg:InvalidArgument", ...
+ "Option 'NewtonTol' will be ignored by this solver");
+ endif
+ if (! isequal (vodeoptions.MaxNewtonIterations,...
vodetemp.MaxNewtonIterations))
- warning ('OdePkg:InvalidArgument', ...
- 'Option "MaxNewtonIterations" will be ignored by this solver');
- end
+ warning ("OdePkg:InvalidArgument", ...
+ "Option 'MaxNewtonIterations' will be ignored by this solver");
+ endif
- %# Implementation of the option Mass has been finished. This option
- %# can be set by the user to another value than default value.
- if (~isempty (vodeoptions.Mass) && isnumeric (vodeoptions.Mass))
- vhavemasshandle = false; vmass = vodeoptions.Mass; %# constant mass
- elseif (isa (vodeoptions.Mass, 'function_handle'))
- vhavemasshandle = true; %# mass defined by a function handle
- else %# no mass matrix - creating a diag-matrix of ones for mass
- vhavemasshandle = false; %# vmass = diag (ones (length (vinit), 1), 0);
- end
+ ## Implementation of the option Mass has been finished. This option
+ ## can be set by the user to another value than default value.
+ if (! isempty (vodeoptions.Mass) && isnumeric (vodeoptions.Mass))
+ vhavemasshandle = false; vmass = vodeoptions.Mass; # constant mass
+ elseif (isa (vodeoptions.Mass, "function_handle"))
+ vhavemasshandle = true; # mass defined by a function handle
+ else # no mass matrix - creating a diag-matrix of ones for mass
+ vhavemasshandle = false; # vmass = diag (ones (length (vinit), 1), 0);
+ endif
- %# Implementation of the option MStateDependence has been finished.
- %# This option can be set by the user to another value than default
- %# value.
- if (strcmp (vodeoptions.MStateDependence, 'none'))
+ ## Implementation of the option MStateDependence has been finished.
+ ## This option can be set by the user to another value than default
+ ## value.
+ if (strcmp (vodeoptions.MStateDependence, "none"))
vmassdependence = false;
else vmassdependence = true;
- end
+ endif
- %# Other options that are not used by this solver. Print a warning
- %# message to tell the user that the option(s) is/are ignored.
- if (~isequal (vodeoptions.MvPattern, vodetemp.MvPattern))
- warning ('OdePkg:InvalidArgument', ...
- 'Option "MvPattern" will be ignored by this solver');
- end
- if (~isequal (vodeoptions.MassSingular, vodetemp.MassSingular))
- warning ('OdePkg:InvalidArgument', ...
- 'Option "MassSingular" will be ignored by this solver');
- end
- if (~isequal (vodeoptions.InitialSlope, vodetemp.InitialSlope))
- warning ('OdePkg:InvalidArgument', ...
- 'Option "InitialSlope" will be ignored by this solver');
- end
- if (~isequal (vodeoptions.MaxOrder, vodetemp.MaxOrder))
- warning ('OdePkg:InvalidArgument', ...
- 'Option "MaxOrder" will be ignored by this solver');
- end
- if (~isequal (vodeoptions.BDF, vodetemp.BDF))
- warning ('OdePkg:InvalidArgument', ...
- 'Option "BDF" will be ignored by this solver');
- end
+ ## Other options that are not used by this solver. Print a warning
+ ## message to tell the user that the option(s) is/are ignored.
+ if (! isequal (vodeoptions.MvPattern, vodetemp.MvPattern))
+ warning ("OdePkg:InvalidArgument", ...
+ "Option 'MvPattern' will be ignored by this solver");
+ endif
+ if (! isequal (vodeoptions.MassSingular, vodetemp.MassSingular))
+ warning ("OdePkg:InvalidArgument", ...
+ "Option 'MassSingular' will be ignored by this solver");
+ endif
+ if (! isequal (vodeoptions.InitialSlope, vodetemp.InitialSlope))
+ warning ("OdePkg:InvalidArgument", ...
+ "Option 'InitialSlope' will be ignored by this solver");
+ endif
+ if (! isequal (vodeoptions.MaxOrder, vodetemp.MaxOrder))
+ warning ("OdePkg:InvalidArgument", ...
+ "Option 'MaxOrder' will be ignored by this solver");
+ endif
+ if (! isequal (vodeoptions.BDF, vodetemp.BDF))
+ warning ("OdePkg:InvalidArgument", ...
+ "Option 'BDF' will be ignored by this solver");
+ endif
- %# Starting the initialisation of the core solver ode23
- vtimestamp = vslot(1,1); %# timestamp = start time
- vtimelength = length (vslot); %# length needed if fixed steps
- vtimestop = vslot(1,vtimelength); %# stop time = last value
- %# 20110611, reported by Nils Strunk
- %# Make it possible to solve equations from negativ to zero,
- %# eg. vres = ode23 (@(t,y) y, [-2 0], 2);
- vdirection = sign (vtimestop - vtimestamp); %# Direction flag
+ ## Starting the initialisation of the core solver ode23
+ vtimestamp = vslot(1, 1); # timestamp = start time
+ vtimelength = length (vslot); # length needed if fixed steps
+ vtimestop = vslot(1, vtimelength); # stop time = last value
+ ## 20110611, reported by Nils Strunk
+ ## Make it possible to solve equations from negativ to zero,
+ ## eg. vres = ode23 (@(t, y) y, [-2 0], 2);
+ vdirection = sign (vtimestop - vtimestamp); # Direction flag
- if (~vstepsizefixed)
+ if (! vstepsizefixed)
if (sign (vodeoptions.InitialStep) == vdirection)
vstepsize = vodeoptions.InitialStep;
- else %# Fix wrong direction of InitialStep.
+ else # Fix wrong direction of InitialStep.
vstepsize = - vodeoptions.InitialStep;
- end
- vminstepsize = (vtimestop - vtimestamp) / (1/eps);
- else %# If step size is given then use the fixed time steps
- vstepsize = vslot(1,2) - vslot(1,1);
+ endif
+ vminstepsize = (vtimestop - vtimestamp) / (1 / eps);
+ else # If step size is given then use the fixed time steps
+ vstepsize = vslot(1, 2) - vslot(1, 1);
vminstepsize = sign (vstepsize) * eps;
- end
+ endif
- vretvaltime = vtimestamp; %# first timestamp output
- vretvalresult = vinit; %# first solution output
+ vretvaltime = vtimestamp; # first timestamp output
+ vretvalresult = vinit; # first solution output
- %# Initialize the OutputFcn
+ ## Initialize the OutputFcn
if (vhaveoutputfunction)
if (vhaveoutputselection) vretout = vretvalresult(vodeoptions.OutputSel);
- else vretout = vretvalresult; end
+ else vretout = vretvalresult;
+ endif
feval (vodeoptions.OutputFcn, vslot.', ...
- vretout.', 'init', vfunarguments{:});
- end
+ vretout.', "init", vfunarguments{:});
+ endif
- %# Initialize the EventFcn
+ ## Initialize the EventFcn
if (vhaveeventfunction)
odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
- vretvalresult.', 'init', vfunarguments{:});
- end
+ vretvalresult.', "init", vfunarguments{:});
+ endif
- vpow = 1/3; %# 20071016, reported by Luis Randez
- va = [ 0, 0, 0; %# The Runge-Kutta-Fehlberg 2(3) coefficients
- 1/2, 0, 0; %# Coefficients proved on 20060827
- -1, 2, 0]; %# See p.91 in Ascher & Petzold
- vb2 = [0; 1; 0]; %# 2nd and 3rd order
- vb3 = [1/6; 2/3; 1/6]; %# b-coefficients
+ vpow = 1 / 3; # 20071016, reported by Luis Randez
+ va = [ 0, 0, 0; # The Runge-Kutta-Fehlberg 2(3) coefficients
+ 1 / 2, 0, 0; # Coefficients proved on 20060827
+ -1, 2, 0]; # See p.91 in Ascher & Petzold
+ vb2 = [0; 1; 0]; # 2nd and 3rd order
+ vb3 = [1 / 6; 2 / 3; 1 / 6]; # b-coefficients
vc = sum (va, 2);
- %# The solver main loop - stop if the endpoint has been reached
- vcntloop = 2; vcntcycles = 1; vu = vinit; vk = vu.' * zeros(1,3);
+ ## The solver main loop - stop if the endpoint has been reached
+ vcntloop = 2; vcntcycles = 1; vu = vinit; vk = vu.' * zeros(1, 3);
vcntiter = 0; vunhandledtermination = true; vcntsave = 2;
while ((vdirection * (vtimestamp) < vdirection * (vtimestop)) && ...
(vdirection * (vstepsize) >= vdirection * (vminstepsize)))
- %# Hit the endpoint of the time slot exactely
+ ## Hit the endpoint of the time slot exactely
if (vdirection * (vtimestamp + vstepsize) > vdirection * vtimestop)
- %# vstepsize = vtimestop - vdirection * vtimestamp;
- %# 20110611, reported by Nils Strunk
- %# The endpoint of the time slot must be hit exactly,
- %# eg. vsol = ode23 (@(t,y) y, [0 -1], 1);
+ ## vstepsize = vtimestop - vdirection * vtimestamp;
+ ## 20110611, reported by Nils Strunk
+ ## The endpoint of the time slot must be hit exactly,
+ ## eg. vsol = ode23 (@(t, y) y, [0 -1], 1);
vstepsize = vdirection * abs (abs (vtimestop) - abs (vtimestamp));
- end
+ endif
- %# Estimate the three results when using this solver
+ ## Estimate the three results when using this solver
for j = 1:3
- vthetime = vtimestamp + vc(j,1) * vstepsize;
- vtheinput = vu.' + vstepsize * vk(:,1:j-1) * va(j,1:j-1).';
- if (vhavemasshandle) %# Handle only the dynamic mass matrix,
- if (vmassdependence) %# constant mass matrices have already
- vmass = feval ... %# been set before (if any)
+ vthetime = vtimestamp + vc(j, 1) * vstepsize;
+ vtheinput = vu.' + vstepsize * vk(:, 1:j - 1) * va(j, 1:j - 1).';
+ if (vhavemasshandle) # Handle only the dynamic mass matrix,
+ if (vmassdependence) # constant mass matrices have already
+ vmass = feval ... # been set before (if any)
(vodeoptions.Mass, vthetime, vtheinput, vfunarguments{:});
- else %# if (vmassdependence == false)
- vmass = feval ... %# then we only have the time argument
+ else # if (vmassdependence == false)
+ vmass = feval ... # then we only have the time argument
(vodeoptions.Mass, vthetime, vfunarguments{:});
- end
- vk(:,j) = vmass \ feval ...
+ endif
+ vk(:, j) = vmass \ feval ...
(vfun, vthetime, vtheinput, vfunarguments{:});
else
- vk(:,j) = feval ...
+ vk(:, j) = feval ...
(vfun, vthetime, vtheinput, vfunarguments{:});
- end
- end
+ endif
+ endfor
- %# Compute the 2nd and the 3rd order estimation
+ ## Compute the 2nd and the 3rd order estimation
y2 = vu.' + vstepsize * (vk * vb2);
y3 = vu.' + vstepsize * (vk * vb3);
if (vhavenonnegative)
vu(vodeoptions.NonNegative) = abs (vu(vodeoptions.NonNegative));
y2(vodeoptions.NonNegative) = abs (y2(vodeoptions.NonNegative));
y3(vodeoptions.NonNegative) = abs (y3(vodeoptions.NonNegative));
- end
+ endif
if (vhaveoutputfunction && vhaverefine)
vSaveVUForRefine = vu;
- end
+ endif
- %# Calculate the absolute local truncation error and the acceptable error
- if (~vstepsizefixed)
- if (~vnormcontrol)
+ ## Calculate the absolute local truncation error and the acceptable error
+ if (! vstepsizefixed)
+ if (! vnormcontrol)
vdelta = abs (y3 - y2);
vtau = max (vodeoptions.RelTol * abs (vu.'), vodeoptions.AbsTol);
else
vdelta = norm (y3 - y2, Inf);
vtau = max (vodeoptions.RelTol * max (norm (vu.', Inf), 1.0), ...
vodeoptions.AbsTol);
- end
- else %# if (vstepsizefixed == true)
+ endif
+ else # if (vstepsizefixed == true)
vdelta = 1; vtau = 2;
- end
+ endif
- %# If the error is acceptable then update the vretval variables
+ ## If the error is acceptable then update the vretval variables
if (all (vdelta <= vtau))
vtimestamp = vtimestamp + vstepsize;
- vu = y3.'; %# MC2001: the higher order estimation as "local extrapolation"
- %# Save the solution every vodeoptions.OutputSave steps
- if (mod (vcntloop-1,vodeoptions.OutputSave) == 0)
- vretvaltime(vcntsave,:) = vtimestamp;
- vretvalresult(vcntsave,:) = vu;
+ vu = y3.'; # MC2001: the higher order estimation as "local extrapolation"
+ ## Save the solution every vodeoptions.OutputSave steps
+ if (mod (vcntloop - 1, vodeoptions.OutputSave) == 0)
+ vretvaltime(vcntsave, :) = vtimestamp;
+ vretvalresult(vcntsave, :) = vu;
vcntsave = vcntsave + 1;
- end
+ endif
vcntloop = vcntloop + 1; vcntiter = 0;
- %# Call plot only if a valid result has been found, therefore this
- %# code fragment has moved here. Stop integration if plot function
- %# returns false
+ ## Call plot only if a valid result has been found, therefore this
+ ## code fragment has moved here. Stop integration if plot function
+ ## returns false
if (vhaveoutputfunction)
- for vcnt = 0:vodeoptions.Refine %# Approximation between told and t
- if (vhaverefine) %# Do interpolation
+ for vcnt = 0:vodeoptions.Refine # Approximation between told and t
+ if (vhaverefine) # Do interpolation
vapproxtime = (vcnt + 1) * vstepsize / (vodeoptions.Refine + 2);
vapproxvals = vSaveVUForRefine.' + vapproxtime * (vk * vb3);
vapproxtime = (vtimestamp - vstepsize) + vapproxtime;
else
vapproxvals = vu.';
vapproxtime = vtimestamp;
- end
+ endif
if (vhaveoutputselection)
vapproxvals = vapproxvals(vodeoptions.OutputSel);
- end
+ endif
vpltret = feval (vodeoptions.OutputFcn, vapproxtime, ...
vapproxvals, [], vfunarguments{:});
- if vpltret %# Leave refinement loop
+ if vpltret # Leave refinement loop
break;
- end
- end
- if (vpltret) %# Leave main loop
+ endif
+ endfor
+ if (vpltret) # Leave main loop
vunhandledtermination = false;
break;
- end
- end
+ endif
+ endif
- %# Call event only if a valid result has been found, therefore this
- %# code fragment has moved here. Stop integration if veventbreak is
- %# true
+ ## Call event only if a valid result has been found, therefore this
+ ## code fragment has moved here. Stop integration if veventbreak is
+ ## true
if (vhaveeventfunction)
vevent = ...
odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
vu(:), [], vfunarguments{:});
- if (~isempty (vevent{1}) && vevent{1} == 1)
- vretvaltime(vcntloop-1,:) = vevent{3}(end,:);
- vretvalresult(vcntloop-1,:) = vevent{4}(end,:);
+ if (! isempty (vevent{1}) && vevent{1} == 1)
+ vretvaltime(vcntloop - 1, :) = vevent{3}(end, :);
+ vretvalresult(vcntloop - 1, :) = vevent{4}(end, :);
vunhandledtermination = false; break;
- end
- end
- end %# If the error is acceptable ...
+ endif
+ endif
+ endif # If the error is acceptable ...
- %# Update the step size for the next integration step
- if (~vstepsizefixed)
- %# 20080425, reported by Marco Caliari
- %# vdelta cannot be negative (because of the absolute value that
- %# has been introduced) but it could be 0, then replace the zeros
- %# with the maximum value of vdelta
+ ## Update the step size for the next integration step
+ if (! vstepsizefixed)
+ ## 20080425, reported by Marco Caliari
+ ## vdelta cannot be negative (because of the absolute value that
+ ## has been introduced) but it could be 0, then replace the zeros
+ ## with the maximum value of vdelta
vdelta(find (vdelta == 0)) = max (vdelta);
- %# It could happen that max (vdelta) == 0 (ie. that the original
- %# vdelta was 0), in that case we double the previous vstepsize
+ ## It could happen that max (vdelta) == 0 (ie. that the original
+ ## vdelta was 0), in that case we double the previous vstepsize
vdelta(find (vdelta == 0)) = max (vtau) .* (0.4 ^ (1 / vpow));
if (vdirection == 1)
@@ -441,96 +446,96 @@
else
vstepsize = max (- vodeoptions.MaxStep, ...
max (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow));
- end
+ endif
- else %# if (vstepsizefixed)
+ else # if (vstepsizefixed)
if (vcntloop <= vtimelength)
- vstepsize = vslot(vcntloop) - vslot(vcntloop-1);
- else %# Get out of the main integration loop
+ vstepsize = vslot(vcntloop) - vslot(vcntloop - 1);
+ else # Get out of the main integration loop
break;
- end
- end
+ endif
+ endif
- %# Update counters that count the number of iteration cycles
- vcntcycles = vcntcycles + 1; %# Needed for cost statistics
- vcntiter = vcntiter + 1; %# Needed to find iteration problems
+ ## Update counters that count the number of iteration cycles
+ vcntcycles = vcntcycles + 1; # Needed for cost statistics
+ vcntiter = vcntiter + 1; # Needed to find iteration problems
- %# Stop solving because the last 1000 steps no successful valid
- %# value has been found
+ # Stop solving because the last 1000 steps no successful valid
+ # value has been found
if (vcntiter >= 5000)
- error (['Solving has not been successful. The iterative', ...
- ' integration loop exited at time t = %f before endpoint at', ...
- ' tend = %f was reached. This happened because the iterative', ...
- ' integration loop does not find a valid solution at this time', ...
- ' stamp. Try to reduce the value of "InitialStep" and/or', ...
- ' "MaxStep" with the command "odeset".\n'], vtimestamp, vtimestop);
- end
+ error (["Solving has not been successful. The iterative", ...
+ " integration loop exited at time t = %f before endpoint at", ...
+ " tend = %f was reached. This happened because the iterative", ...
+ " integration loop does not find a valid solution at this time", ...
+ " stamp. Try to reduce the value of 'InitialStep' and/or", ...
+ " 'MaxStep' with the command 'odeset'.\n"], vtimestamp, vtimestop);
+ endif
- end %# The main loop
+ endwhile # The main loop
- %# Check if integration of the ode has been successful
+ ## Check if integration of the ode has been successful
if (vdirection * vtimestamp < vdirection * vtimestop)
if (vunhandledtermination == true)
- error ('OdePkg:InvalidArgument', ...
- ['Solving has not been successful. The iterative', ...
- ' integration loop exited at time t = %f', ...
- ' before endpoint at tend = %f was reached. This may', ...
- ' happen if the stepsize grows smaller than defined in', ...
- ' vminstepsize. Try to reduce the value of "InitialStep" and/or', ...
- ' "MaxStep" with the command "odeset".\n'], vtimestamp, vtimestop);
+ error ("OdePkg:InvalidArgument", ...
+ ["Solving has not been successful. The iterative", ...
+ " integration loop exited at time t = %f", ...
+ " before endpoint at tend = %f was reached. This may", ...
+ " happen if the stepsize grows smaller than defined in", ...
+ " vminstepsize. Try to reduce the value of 'InitialStep' and/or", ...
+ " 'MaxStep' with the command 'odeset'.\n"], vtimestamp, vtimestop);
else
- warning ('OdePkg:InvalidArgument', ...
- ['Solver has been stopped by a call of "break" in', ...
- ' the main iteration loop at time t = %f before endpoint at', ...
- ' tend = %f was reached. This may happen because the @odeplot', ...
- ' function returned "true" or the @event function returned "true".'], ...
+ warning ("OdePkg:InvalidArgument", ...
+ ["Solver has been stopped by a call of 'break' in", ...
+ " the main iteration loop at time t = %f before endpoint at", ...
+ " tend = %f was reached. This may happen because the @odeplot", ...
+ " function returned 'true' or the @event function returned 'true'."], ...
vtimestamp, vtimestop);
- end
- end
+ endif
+ endif
- %# Postprocessing, do whatever when terminating integration algorithm
- if (vhaveoutputfunction) %# Cleanup plotter
+ ## Postprocessing, do whatever when terminating integration algorithm
+ if (vhaveoutputfunction) # Cleanup plotter
feval (vodeoptions.OutputFcn, vtimestamp, ...
- vu.', 'done', vfunarguments{:});
- end
- if (vhaveeventfunction) %# Cleanup event function handling
+ vu.', "done", vfunarguments{:});
+ endif
+ if (vhaveeventfunction) # Cleanup event function handling
odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
- vu.', 'done', vfunarguments{:});
- end
- %# Save the last step, if not already saved
- if (mod (vcntloop-2,vodeoptions.OutputSave) ~= 0)
- vretvaltime(vcntsave,:) = vtimestamp;
- vretvalresult(vcntsave,:) = vu;
- end
+ vu.', "done", vfunarguments{:});
+ endif
+ ## Save the last step, if not already saved
+ if (mod (vcntloop - 2, vodeoptions.OutputSave) != 0)
+ vretvaltime(vcntsave, :) = vtimestamp;
+ vretvalresult(vcntsave, :) = vu;
+ endif
- %# Print additional information if option Stats is set
- if (strcmp (vodeoptions.Stats, 'on'))
+ ## Print additional information if option Stats is set
+ if (strcmp (vodeoptions.Stats, "on"))
vhavestats = true;
- vnsteps = vcntloop-2; %# vcntloop from 2..end
- vnfailed = (vcntcycles-1)-(vcntloop-2)+1; %# vcntcycl from 1..end
- vnfevals = 3*(vcntcycles-1); %# number of ode evaluations
- vndecomps = 0; %# number of LU decompositions
- vnpds = 0; %# number of partial derivatives
- vnlinsols = 0; %# no. of solutions of linear systems
- %# Print cost statistics if no output argument is given
+ vnsteps = vcntloop - 2; # vcntloop from 2..end
+ vnfailed = (vcntcycles - 1) - (vcntloop - 2) + 1; # vcntcycl from 1..end
+ vnfevals = 3 * (vcntcycles - 1); # number of ode evaluations
+ vndecomps = 0; # number of LU decompositions
+ vnpds = 0; # number of partial derivatives
+ vnlinsols = 0; # no. of solutions of linear systems
+ ## Print cost statistics if no output argument is given
if (nargout == 0)
- vmsg = fprintf (1, 'Number of successful steps: %d\n', vnsteps);
- vmsg = fprintf (1, 'Number of failed attempts: %d\n', vnfailed);
- vmsg = fprintf (1, 'Number of function calls: %d\n', vnfevals);
- end
+ vmsg = fprintf (1, "Number of successful steps: %d\n", vnsteps);
+ vmsg = fprintf (1, "Number of failed attempts: %d\n", vnfailed);
+ vmsg = fprintf (1, "Number of function calls: %d\n", vnfevals);
+ endif
else
vhavestats = false;
- end
+ endif
- if (nargout == 1) %# Sort output variables, depends on nargout
- varargout{1}.x = vretvaltime; %# Time stamps are saved in field x
- varargout{1}.y = vretvalresult; %# Results are saved in field y
- varargout{1}.solver = 'ode23'; %# Solver name is saved in field solver
+ if (nargout == 1) # Sort output variables, depends on nargout
+ varargout{1}.x = vretvaltime; # Time stamps are saved in field x
+ varargout{1}.y = vretvalresult; # Results are saved in field y
+ varargout{1}.solver = "ode23"; # Solver name is saved in field solver
if (vhaveeventfunction)
- varargout{1}.ie = vevent{2}; %# Index info which event occured
- varargout{1}.xe = vevent{3}; %# Time info when an event occured
- varargout{1}.ye = vevent{4}; %# Results when an event occured
- end
+ varargout{1}.ie = vevent{2}; # Index info which event occured
+ varargout{1}.xe = vevent{3}; # Time info when an event occured
+ varargout{1}.ye = vevent{4}; # Results when an event occured
+ endif
if (vhavestats)
varargout{1}.stats = struct;
varargout{1}.stats.nsteps = vnsteps;
@@ -539,215 +544,215 @@
varargout{1}.stats.npds = vnpds;
varargout{1}.stats.ndecomps = vndecomps;
varargout{1}.stats.nlinsols = vnlinsols;
- end
+ endif
elseif (nargout == 2)
- varargout{1} = vretvaltime; %# Time stamps are first output argument
- varargout{2} = vretvalresult; %# Results are second output argument
+ varargout{1} = vretvaltime; # Time stamps are first output argument
+ varargout{2} = vretvalresult; # Results are second output argument
elseif (nargout == 5)
- varargout{1} = vretvaltime; %# Same as (nargout == 2)
- varargout{2} = vretvalresult; %# Same as (nargout == 2)
- varargout{3} = []; %# LabMat doesn't accept lines like
- varargout{4} = []; %# varargout{3} = varargout{4} = [];
+ varargout{1} = vretvaltime; # Same as (nargout == 2)
+ varargout{2} = vretvalresult; # Same as (nargout == 2)
+ varargout{3} = []; # LabMat doesn't accept lines like
+ varargout{4} = []; # varargout{3} = varargout{4} = [];
varargout{5} = [];
if (vhaveeventfunction)
- varargout{3} = vevent{3}; %# Time info when an event occured
- varargout{4} = vevent{4}; %# Results when an event occured
- varargout{5} = vevent{2}; %# Index info which event occured
- end
- end
-end
+ varargout{3} = vevent{3}; # Time info when an event occured
+ varargout{4} = vevent{4}; # Results when an event occured
+ varargout{5} = vevent{2}; # Index info which event occured
+ endif
+ endif
+endfunction
-%! # We are using the "Van der Pol" implementation for all tests that
-%! # are done for this function. We also define a Jacobian, Events,
-%! # pseudo-Mass implementation. For further tests we also define a
-%! # reference solution (computed at high accuracy) and an OutputFcn
-%!function [ydot] = fpol (vt, vy, varargin) %# The Van der Pol
-%! ydot = [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
-%!function [vjac] = fjac (vt, vy, varargin) %# its Jacobian
-%! vjac = [0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2];
-%!function [vjac] = fjcc (vt, vy, varargin) %# sparse type
-%! vjac = sparse ([0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2]);
+%!# We are using the "Van der Pol" implementation for all tests that
+%!# are done for this function. We also define a Jacobian, Events,
+%!# pseudo-Mass implementation. For further tests we also define a
+%!# reference solution (computed at high accuracy) and an OutputFcn
+%!function [ydot] = fpol (vt, vy, varargin) # The Van der Pol
+%! ydot = [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
+%!function [vjac] = fjac (vt, vy, varargin) # its Jacobian
+%! vjac = [0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2];
+%!function [vjac] = fjcc (vt, vy, varargin) # sparse type
+%! vjac = sparse ([0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2]);
%!function [vval, vtrm, vdir] = feve (vt, vy, varargin)
-%! vval = fpol (vt, vy, varargin); %# We use the derivatives
-%! vtrm = zeros (2,1); %# that's why component 2
-%! vdir = ones (2,1); %# seems to not be exact
+%! vval = fpol (vt, vy, varargin); # We use the derivatives
+%! vtrm = zeros (2,1); # that's why component 2
+%! vdir = ones (2,1); # seems to not be exact
%!function [vval, vtrm, vdir] = fevn (vt, vy, varargin)
-%! vval = fpol (vt, vy, varargin); %# We use the derivatives
-%! vtrm = ones (2,1); %# that's why component 2
-%! vdir = ones (2,1); %# seems to not be exact
+%! vval = fpol (vt, vy, varargin); # We use the derivatives
+%! vtrm = ones (2, 1); # that's why component 2
+%! vdir = ones (2, 1); # seems to not be exact
%!function [vmas] = fmas (vt, vy)
-%! vmas = [1, 0; 0, 1]; %# Dummy mass matrix for tests
+%! vmas = [1, 0; 0, 1]; # Dummy mass matrix for tests
%!function [vmas] = fmsa (vt, vy)
-%! vmas = sparse ([1, 0; 0, 1]); %# A sparse dummy matrix
-%!function [vref] = fref () %# The computed reference sol
-%! vref = [0.32331666704577, -1.83297456798624];
+%! vmas = sparse ([1, 0; 0, 1]); # A sparse dummy matrix
+%!function [vref] = fref () # The computed reference sol
+%! vref = [0.32331666704577, -1.83297456798624];
%!function [vout] = fout (vt, vy, vflag, varargin)
-%! if (regexp (char (vflag), 'init') == 1)
-%! if (any (size (vt) ~= [2, 1])) error ('"fout" step "init"'); end
-%! elseif (isempty (vflag))
-%! if (any (size (vt) ~= [1, 1])) error ('"fout" step "calc"'); end
-%! vout = false;
-%! elseif (regexp (char (vflag), 'done') == 1)
-%! if (any (size (vt) ~= [1, 1])) error ('"fout" step "done"'); end
-%! else error ('"fout" invalid vflag');
-%! end
+%! if (regexp (char (vflag), "init") == 1)
+%! if (any (size (vt) != [2, 1])) error ("'fout' step 'init'"); end
+%! elseif (isempty (vflag))
+%! if (any (size (vt) != [1, 1])) error ("'fout' step 'calc'"); end
+%! vout = false;
+%! elseif (regexp (char (vflag), "done") == 1)
+%! if (any (size (vt) != [1, 1])) error ("'fout' step 'done'"); end
+%! else error ("'fout' invalid vflag");
+%! end
%!
-%! %# Turn off output of warning messages for all tests, turn them on
-%! %# again if the last test is called
-%!error %# input argument number one
-%! warning ('off', 'OdePkg:InvalidArgument');
-%! B = ode23 (1, [0 25], [3 15 1]);
-%!error %# input argument number two
-%! B = ode23 (@fpol, 1, [3 15 1]);
-%!error %# input argument number three
-%! B = ode23 (@flor, [0 25], 1);
-%!test %# one output argument
-%! vsol = ode23 (@fpol, [0 2], [2 0]);
-%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
-%! assert (isfield (vsol, 'solver'));
-%! assert (vsol.solver, 'ode23');
-%!test %# two output arguments
-%! [vt, vy] = ode23 (@fpol, [0 2], [2 0]);
-%! assert ([vt(end), vy(end,:)], [2, fref], 1e-3);
-%!test %# five output arguments and no Events
-%! [vt, vy, vxe, vye, vie] = ode23 (@fpol, [0 2], [2 0]);
-%! assert ([vt(end), vy(end,:)], [2, fref], 1e-3);
-%! assert ([vie, vxe, vye], []);
-%!test %# anonymous function instead of real function
-%! fvdb = @(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
-%! vsol = ode23 (fvdb, [0 2], [2 0]);
-%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
-%!test %# extra input arguments passed through
-%! vsol = ode23 (@fpol, [0 2], [2 0], 12, 13, 'KL');
-%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
-%!test %# empty OdePkg structure *but* extra input arguments
-%! vopt = odeset;
-%! vsol = ode23 (@fpol, [0 2], [2 0], vopt, 12, 13, 'KL');
-%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
-%!error %# strange OdePkg structure
-%! vopt = struct ('foo', 1);
-%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
-%!test %# Solve vdp in fixed step sizes
-%! vsol = ode23 (@fpol, [0:0.1:2], [2 0]);
-%! assert (vsol.x(:), [0:0.1:2]');
-%! assert (vsol.y(end,:), fref, 1e-3);
-%!test %# Solve in backward direction starting at t=0
-%! vref = [-1.205364552835178, 0.951542399860817];
-%! vsol = ode23 (@fpol, [0 -2], [2 0]);
-%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-3);
-%!test %# Solve in backward direction starting at t=2
-%! vref = [-1.205364552835178, 0.951542399860817];
-%! vsol = ode23 (@fpol, [2 -2], fref);
-%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-3);
-%!test %# Solve another anonymous function in backward direction
-%! vref = [-1, 0.367879437558975];
-%! vsol = ode23 (@(t,y) y, [0 -1], 1);
-%! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3);
-%!test %# Solve another anonymous function below zero
-%! vref = [0, 14.77810590694212];
-%! vsol = ode23 (@(t,y) y, [-2 0], 2);
-%! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3);
-%!test %# AbsTol option
-%! vopt = odeset ('AbsTol', 1e-5);
-%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
-%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
-%!test %# AbsTol and RelTol option
-%! vopt = odeset ('AbsTol', 1e-8, 'RelTol', 1e-8);
-%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
-%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
-%!test %# RelTol and NormControl option -- higher accuracy
-%! vopt = odeset ('RelTol', 1e-8, 'NormControl', 'on');
-%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
-%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-6);
-%!test %# Keeps initial values while integrating
-%! vopt = odeset ('NonNegative', 2);
-%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
-%! assert ([vsol.x(end), vsol.y(end,:)], [2, 2, 0], 1e-1);
-%!test %# Details of OutputSel and Refine can't be tested
-%! vopt = odeset ('OutputFcn', @fout, 'OutputSel', 1, 'Refine', 5);
-%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
-%!test %# Details of OutputSave can't be tested
-%! vopt = odeset ('OutputSave', 1, 'OutputSel', 1);
-%! vsla = ode23 (@fpol, [0 2], [2 0], vopt);
-%! vopt = odeset ('OutputSave', 2);
-%! vslb = ode23 (@fpol, [0 2], [2 0], vopt);
-%! assert (length (vsla.x) > length (vslb.x))
-%!test %# Stats must add further elements in vsol
-%! vopt = odeset ('Stats', 'on');
-%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
-%! assert (isfield (vsol, 'stats'));
-%! assert (isfield (vsol.stats, 'nsteps'));
-%!test %# InitialStep option
-%! vopt = odeset ('InitialStep', 1e-8);
-%! vsol = ode23 (@fpol, [0 0.2], [2 0], vopt);
-%! assert ([vsol.x(2)-vsol.x(1)], [1e-8], 1e-9);
-%!test %# MaxStep option
-%! vopt = odeset ('MaxStep', 1e-2);
-%! vsol = ode23 (@fpol, [0 0.2], [2 0], vopt);
-%! assert ([vsol.x(5)-vsol.x(4)], [1e-2], 1e-2);
-%!test %# Events option add further elements in vsol
-%! vopt = odeset ('Events', @feve);
-%! vsol = ode23 (@fpol, [0 10], [2 0], vopt);
-%! assert (isfield (vsol, 'ie'));
-%! assert (vsol.ie, [2; 1; 2; 1]);
-%! assert (isfield (vsol, 'xe'));
-%! assert (isfield (vsol, 'ye'));
-%!test %# Events option, now stop integration
-%! vopt = odeset ('Events', @fevn, 'NormControl', 'on');
-%! vsol = ode23 (@fpol, [0 10], [2 0], vopt);
-%! assert ([vsol.ie, vsol.xe, vsol.ye], ...
-%! [2.0, 2.496110, -0.830550, -2.677589], 1e-3);
-%!test %# Events option, five output arguments
-%! vopt = odeset ('Events', @fevn, 'NormControl', 'on');
-%! [vt, vy, vxe, vye, vie] = ode23 (@fpol, [0 10], [2 0], vopt);
-%! assert ([vie, vxe, vye], ...
-%! [2.0, 2.496110, -0.830550, -2.677589], 1e-3);
-%!test %# Jacobian option
-%! vopt = odeset ('Jacobian', @fjac);
-%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
-%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
-%!test %# Jacobian option and sparse return value
-%! vopt = odeset ('Jacobian', @fjcc);
-%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
-%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
+%!## Turn off output of warning messages for all tests, turn them on
+%!## again if the last test is called
+%!error # input argument number one
+%! warning ("off", "OdePkg:InvalidArgument");
+%! B = ode23 (1, [0 25], [3 15 1]);
+%!error # input argument number two
+%! B = ode23 (@fpol, 1, [3 15 1]);
+%!error # input argument number three
+%! B = ode23 (@flor, [0 25], 1);
+%!test # one output argument
+%! vsol = ode23 (@fpol, [0 2], [2 0]);
+%! assert ([vsol.x(end), vsol.y(end, :)], [2, fref], 1e-3);
+%! assert (isfield (vsol, "solver"));
+%! assert (vsol.solver, "ode23");
+%!test # two output arguments
+%! [vt, vy] = ode23 (@fpol, [0 2], [2 0]);
+%! assert ([vt(end), vy(end, :)], [2, fref], 1e-3);
+%!test # five output arguments and no Events
+%! [vt, vy, vxe, vye, vie] = ode23 (@fpol, [0 2], [2 0]);
+%! assert ([vt(end), vy(end, :)], [2, fref], 1e-3);
+%! assert ([vie, vxe, vye], []);
+%!test # anonymous function instead of real function
+%! fvdb = @(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
+%! vsol = ode23 (fvdb, [0 2], [2 0]);
+%! assert ([vsol.x(end), vsol.y(end, :)], [2, fref], 1e-3);
+%!test # extra input arguments passed through
+%! vsol = ode23 (@fpol, [0 2], [2 0], 12, 13, "KL");
+%! assert ([vsol.x(end), vsol.y(end, :)], [2, fref], 1e-3);
+%!test # empty OdePkg structure *but* extra input arguments
+%! vopt = odeset;
+%! vsol = ode23 (@fpol, [0 2], [2 0], vopt, 12, 13, "KL");
+%! assert ([vsol.x(end), vsol.y(end, :)], [2, fref], 1e-3);
+%!error # strange OdePkg structure
+%! vopt = struct ("foo", 1);
+%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
+%!test # Solve vdp in fixed step sizes
+%! vsol = ode23 (@fpol, [0:0.1:2], [2 0]);
+%! assert (vsol.x(:), [0:0.1:2]');
+%! assert (vsol.y(end, :), fref, 1e-3);
+%!test # Solve in backward direction starting at t=0
+%! vref = [-1.205364552835178, 0.951542399860817];
+%! vsol = ode23 (@fpol, [0 -2], [2 0]);
+%! assert ([vsol.x(end), vsol.y(end, :)], [-2, vref], 1e-3);
+%!test # Solve in backward direction starting at t=2
+%! vref = [-1.205364552835178, 0.951542399860817];
+%! vsol = ode23 (@fpol, [2 -2], fref);
+%! assert ([vsol.x(end), vsol.y(end, :)], [-2, vref], 1e-3);
+%!test # Solve another anonymous function in backward direction
+%! vref = [-1, 0.367879437558975];
+%! vsol = ode23 (@(t,y) y, [0 -1], 1);
+%! assert ([vsol.x(end), vsol.y(end, :)], vref, 1e-3);
+%!test # Solve another anonymous function below zero
+%! vref = [0, 14.77810590694212];
+%! vsol = ode23 (@(t,y) y, [-2 0], 2);
+%! assert ([vsol.x(end), vsol.y(end, :)], vref, 1e-3);
+%!test # AbsTol option
+%! vopt = odeset ("AbsTol", 1e-5);
+%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
+%! assert ([vsol.x(end), vsol.y(end, :)], [2, fref], 1e-3);
+%!test # AbsTol and RelTol option
+%! vopt = odeset ("AbsTol", 1e-8, "RelTol", 1e-8);
+%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
+%! assert ([vsol.x(end), vsol.y(end, :)], [2, fref], 1e-3);
+%!test # RelTol and NormControl option -- higher accuracy
+%! vopt = odeset ("RelTol", 1e-8, "NormControl", "on");
+%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
+%! assert ([vsol.x(end), vsol.y(end, :)], [2, fref], 1e-6);
+%!test # Keeps initial values while integrating
+%! vopt = odeset ("NonNegative", 2);
+%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
+%! assert ([vsol.x(end), vsol.y(end, :)], [2, 2, 0], 1e-1);
+%!test # Details of OutputSel and Refine can't be tested
+%! vopt = odeset ("OutputFcn", @fout, "OutputSel", 1, "Refine", 5);
+%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
+%!test # Details of OutputSave can't be tested
+%! vopt = odeset ("OutputSave", 1, "OutputSel", 1);
+%! vsla = ode23 (@fpol, [0 2], [2 0], vopt);
+%! vopt = odeset ("OutputSave", 2);
+%! vslb = ode23 (@fpol, [0 2], [2 0], vopt);
+%! assert (length (vsla.x) > length (vslb.x))
+%!test # Stats must add further elements in vsol
+%! vopt = odeset ("Stats", "on");
+%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
+%! assert (isfield (vsol, "stats"));
+%! assert (isfield (vsol.stats, "nsteps"));
+%!test # InitialStep option
+%! vopt = odeset ("InitialStep", 1e-8);
+%! vsol = ode23 (@fpol, [0 0.2], [2 0], vopt);
+%! assert ([vsol.x(2)-vsol.x(1)], [1e-8], 1e-9);
+%!test # MaxStep option
+%! vopt = odeset ("MaxStep", 1e-2);
+%! vsol = ode23 (@fpol, [0 0.2], [2 0], vopt);
+%! assert ([vsol.x(5)-vsol.x(4)], [1e-2], 1e-2);
+%!test # Events option add further elements in vsol
+%! vopt = odeset ("Events", @feve);
+%! vsol = ode23 (@fpol, [0 10], [2 0], vopt);
+%! assert (isfield (vsol, "ie"));
+%! assert (vsol.ie, [2; 1; 2; 1]);
+%! assert (isfield (vsol, "xe"));
+%! assert (isfield (vsol, "ye"));
+%!test # Events option, now stop integration
+%! vopt = odeset ("Events", @fevn, "NormControl", "on");
+%! vsol = ode23 (@fpol, [0 10], [2 0], vopt);
+%! assert ([vsol.ie, vsol.xe, vsol.ye], ...
+%! [2.0, 2.496110, -0.830550, -2.677589], 1e-3);
+%!test # Events option, five output arguments
+%! vopt = odeset ("Events", @fevn, "NormControl", "on");
+%! [vt, vy, vxe, vye, vie] = ode23 (@fpol, [0 10], [2 0], vopt);
+%! assert ([vie, vxe, vye], ...
+%! [2.0, 2.496110, -0.830550, -2.677589], 1e-3);
+%!test # Jacobian option
+%! vopt = odeset ("Jacobian", @fjac);
+%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
+%! assert ([vsol.x(end), vsol.y(end, :)], [2, fref], 1e-3);
+%!test # Jacobian option and sparse return value
+%! vopt = odeset ("Jacobian", @fjcc);
+%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
+%! assert ([vsol.x(end), vsol.y(end, :)], [2, fref], 1e-3);
%!
-%! %# test for JPattern option is missing
-%! %# test for Vectorized option is missing
-%! %# test for NewtonTol option is missing
-%! %# test for MaxNewtonIterations option is missing
+%!## test for JPattern option is missing
+%!## test for Vectorized option is missing
+%!## test for NewtonTol option is missing
+%!## test for MaxNewtonIterations option is missing
%!
-%!test %# Mass option as function
-%! vopt = odeset ('Mass', @fmas);
-%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
-%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
-%!test %# Mass option as matrix
-%! vopt = odeset ('Mass', eye (2,2));
-%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
-%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
-%!test %# Mass option as sparse matrix
-%! vopt = odeset ('Mass', sparse (eye (2,2)));
-%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
-%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
-%!test %# Mass option as function and sparse matrix
-%! vopt = odeset ('Mass', @fmsa);
-%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
-%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
-%!test %# Mass option as function and MStateDependence
-%! vopt = odeset ('Mass', @fmas, 'MStateDependence', 'strong');
-%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
-%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
-%!test %# Set BDF option to something else than default
-%! vopt = odeset ('BDF', 'on');
-%! [vt, vy] = ode23 (@fpol, [0 2], [2 0], vopt);
-%! assert ([vt(end), vy(end,:)], [2, fref], 1e-3);
+%!test # Mass option as function
+%! vopt = odeset ("Mass", @fmas);
+%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
+%! assert ([vsol.x(end), vsol.y(end, :)], [2, fref], 1e-3);
+%!test # Mass option as matrix
+%! vopt = odeset ("Mass", eye (2,2));
+%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
+%! assert ([vsol.x(end), vsol.y(end, :)], [2, fref], 1e-3);
+%!test # Mass option as sparse matrix
+%! vopt = odeset ("Mass", sparse (eye (2,2)));
+%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
+%! assert ([vsol.x(end), vsol.y(end, :)], [2, fref], 1e-3);
+%!test # Mass option as function and sparse matrix
+%! vopt = odeset ("Mass", @fmsa);
+%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
+%! assert ([vsol.x(end), vsol.y(end, :)], [2, fref], 1e-3);
+%!test # Mass option as function and MStateDependence
+%! vopt = odeset ("Mass", @fmas, "MStateDependence", "strong");
+%! vsol = ode23 (@fpol, [0 2], [2 0], vopt);
+%! assert ([vsol.x(end), vsol.y(end, :)], [2, fref], 1e-3);
+%!test # Set BDF option to something else than default
+%! vopt = odeset ("BDF", "on");
+%! [vt, vy] = ode23 (@fpol, [0 2], [2 0], vopt);
+%! assert ([vt(end), vy(end, :)], [2, fref], 1e-3);
%!
-%! %# test for MvPattern option is missing
-%! %# test for InitialSlope option is missing
-%! %# test for MaxOrder option is missing
+%!## test for MvPattern option is missing
+%!## test for InitialSlope option is missing
+%!## test for MaxOrder option is missing
%!
-%! warning ('on', 'OdePkg:InvalidArgument');
+%! warning ("on", "OdePkg:InvalidArgument");
-%# Local Variables: ***
-%# mode: octave ***
-%# End: ***
+## Local Variables: ***
+## mode: octave ***
+## End: ***