--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/ode/decic.m	Tue Aug 23 02:31:51 2016 +0200
@@ -0,0 +1,223 @@
+## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it>
+##
+## This file is part of Octave.
+##
+## Octave 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 3 of the License, or (at
+## your option) any later version.
+##
+## Octave 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 Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {} {[@var{y0_new}, @var{yp0_new}] =} decic (@var{fun}, @var{t0}, @var{y0}, @var{fixed_y0}, @var{yp0}, @var{fixed_yp0})
+## @deftypefnx  {} {[@var{y0_new}, @var{yp0_new}] =} decic (@var{fun}, @var{t0}, @var{y0}, @var{fixed_y0}, @var{yp0}, @var{fixed_yp0}, @var{options})
+## @deftypefnx  {} {[@var{y0_new}, @var{yp0_new}, @var{resnorm}] =} decic (@dots{})
+##
+## Compute consistent initial conditions @var{y0_new} and @var{yp0_new}, given
+## initial guesses @var{y0} and @var{yp0}. Choose a maximum of
+## @code{length(@var{y0})} components between @var{fixed_y0} and
+## @var{fixed_yp0} as fixed values.
+##
+## @var{fun} is a function handle. The function must accept three inputs where
+## the first is time @var{t}, the second is a column vector of unknowns @var{y}
+## and the third is a column vector of unknowns @var{yp}.
+##
+## @var{t0} is the initial time such that @code{@var{fun}(@var{t0},
+## @var{y0_new}, @var{yp0_new}) = 0}, specified as a scalar. 
+##
+## @var{y0} is a vector used as initial guess for @var{y}.
+##
+## @var{fixed_y0} is a vector which specifies the components of @var{y0} to be
+## hold fixed. Choose a maximum of @code{length(@var{y0})} components between
+## @var{fixed_y0} and @var{fixed_yp0} as fixed values.
+## Set @var{fixed_y0}(i) component to 1 if you want to fix the value of
+## @var{y0}(i).
+## Set @var{fixed_y0}(i) component to 0 if you want to allow the value of
+## @var{y0}(i) to change.
+##
+## @var{yp0} is a vector used as initial guess for @var{yp}.
+##
+## @var{fixed_yp0} is a vector which specifies the components of @var{yp0} to
+## be hold fixed. Choose a maximum of @code{length(@var{yp0})} components
+## between @var{fixed_y0} and @var{fixed_yp0} as fixed values.
+## Set @var{fixed_yp0}(i) component to 1 if you want to fix the value of
+## @var{yp0}(i).
+## Set @var{fixed_yp0}(i) component to 0 if you want to allow the value of
+## @var{yp0}(i) to change.
+##
+## The optional seventh argument @var{options} is a structure array.
+## Use @code{odeset} to generate this structure. The relevant options are
+## @code{RelTol} and @code{AbsTol} which specify the error thresholds used to
+## compute the initial conditions. 
+##
+## The function typically returns two outputs. Variable @var{y0_new} is a
+## column vector and contains the consistent initial value of y.  The
+## output @var{yp0_new} is a column vector and contains the consistent initial
+## value of yp.
+##
+## The optional third output @var{resnorm} is the vector of norm of the 
+## residuals. If @var{resnorm} is small, @code{decic} has successfully computed
+## the initial conditions. If the value of  @var{resnorm} is large, use
+## @code{RelTol} and @code{AbsTol} to adjust it.  
+##
+## Example: Compute initial conditions of @nospell{Robetson}'s equations:
+##
+## @example
+## @group
+## function res = robertsidae(@var{t}, @var{y}, @var{yp})
+## res = [-(@var{yp}(1) + 0.04*@var{y}(1) - 1e4*@var{y}(2)*@var{y}(3));
+##        -(@var{yp}(2) - 0.04*@var{y}(1) + 1e4*@var{y}(2)*@var{y}(3) +
+##        3e7*@var{y}(2)^2);
+##        @var{y}(1) + @var{y}(2) + @var{y}(3) - 1];
+## endfunction
+## @end group
+## [@var{y0_new},@var{yp0_new}] = decic (@@robertsidae, 0, [1; 0; 0], [1; 1; 0],
+## [-1e-4; 1; 0], [0; 0; 0]);
+## @end example
+## @seealso{odeset}
+## @end deftypefn
+
+function [y0_new, yp0_new, resnrm] = decic (odefun, t0, y0, fixed_y0, yp0,
+                                            fixed_yp0, options)
+
+  if (nargin < 6 || nargin > 7 || nargout > 3)
+    print_usage ();
+  endif
+
+  #Check input
+  if (! isa (odefun, "function_handle"))
+    error ("Octave:invalid-input-arg",
+           "decic: FUN must be a valid function handle");
+  endif
+
+  if (! isnumeric (t0) || numel (t0) != 1)
+    error ("Octave:invalid-input-arg",
+           "decic: INIT must be a numeric scalar value");
+  endif
+
+  if (! isnumeric (y0) || ! isvector (y0) || ! isnumeric (fixed_y0) ||
+      ! isvector (fixed_y0) || ! isnumeric (yp0) || ! isvector (yp0)|| 
+      ! isnumeric (fixed_yp0) || ! isvector (fixed_yp0))
+   error ("Octave:invalid-input-arg",
+           "decic: y0, fixed_y0, yp0 and fixed_yp0 must be numeric vectors");
+  
+  elseif (! isequal (numel (y0), numel (fixed_y0), numel (yp0),
+          numel (fixed_yp0)))
+    error ("Octave:invalid-input-arg",
+           "decic: length of y0, fixed_y0, yp0 and fixed_yp0 must be equal");
+  endif
+
+  for i = 1:numel (y0)
+    if (! (fixed_y0(i) == 0 || fixed_y0(i) == 1) || ! (fixed_yp0(i) == 0
+        || fixed_yp0(i) == 1))
+      error ("Octave:invalid-input-arg",
+             "decic: fixed_y0 and fixed_yp0 must be boolean vectors");
+    endif
+  endfor
+  
+  n  = numel (y0);
+  nl = sum (~fixed_y0);
+  nu = sum (~fixed_yp0);
+ 
+  if (n - nl - nu > 0)
+    error ("Octave:invalid-input-arg",
+           "decic: you cannot fix more than length(y0) components");
+  endif
+
+  #Set default value
+  TolFun = 0;
+  TolX   = eps;
+  
+  #Check AbsTol and RelTol
+  if (nargin == 7)
+    if (! isempty (options.AbsTol))
+      if (! isscalar (options.AbsTol))
+        error ("Octave:invalid-input-arg",
+               "decic: AbsTol must be a scalar value");
+      else
+        TolFun = options.AbsTol; 
+      endif
+    endif
+
+    if (! isempty (options.RelTol))
+      if (! isscalar (options.RelTol))
+        error ("Octave:invalid-input-arg",
+               "decic: RelTol must be a scalar value");
+      else
+        TolX = options.RelTol; 
+      endif
+    endif
+  endif
+
+  x0 = [y0(~fixed_y0); yp0(~fixed_yp0)];
+  opt = optimset ('tolfun', TolFun, 'tolx', TolX, 'display', 'iter-detailed');
+  x = fminunc (@(x) objective (x, t0, y0, fixed_y0, yp0, fixed_yp0, nl, nu,
+                               odefun),
+               x0, opt);
+  
+  y0_new  = y0;
+  yp0_new = yp0;
+  
+  y0_new(~fixed_y0)   = x(1:nl);
+  yp0_new(~fixed_yp0) = x(nl+1:nl+nu);
+  resnrm = odefun (t0, y0_new, yp0_new);
+ 
+endfunction
+
+function res = objective (x, t0, y0, fixed_y0, yp0,
+                          fixed_yp0, nl, nu, odefun)
+  y = y0;
+  y(~fixed_y0) = x(1:nl);
+  yp = yp0;
+  yp(~fixed_yp0) = x(nl+1:nl+nu);
+  res = sqrt (sum (odefun (t0, y, yp) .^ 2));
+endfunction
+
+%!function res = rob (t, y, yp)
+%! res =[-(yp(1) + 0.04*y(1) - 1e4*y(2)*y(3)); -(yp(2) - 0.04*y(1) +
+%!      + 1e4*y(2)*y(3) + 3e7*y(2)^2); y(1) + y(2) + y(3) - 1];
+%!endfunction
+
+%!test  # Without options
+%! ref1 = [1;0;0]; 
+%! ref2 = [-4e-2; 4e-2; 0];
+%! [ynew,ypnew] = decic (@rob,0,[1;0;0],[1;1;0],[23;110;0],[0;0;1]);
+%! assert ([ynew(1:end), ypnew(1:end)], [ref1(1:end), ref2(1:end)], 1e-10);
+%!test  # With options
+%! ref1 = [1;0;0]; 
+%! ref2 = [-4e-2; 4e-2; 0];
+%! opt = odeset ("AbsTol", 1e-8, "RelTol", 1e-4);
+%! [ynew,ypnew] = decic (@rob,0,[1;0;0],[1;1;0],[23;110;0],[0;0;1],opt);
+%! assert ([ynew(1:end), ypnew(1:end)], [ref1(1:end), ref2(1:end)], 1e-5);
+
+## Test input validation
+%!error decic ()
+%!error decic (1)
+%!error decic (1,2)
+%!error decic (1,2,3,4,5,6,7,8)
+%!error <FUN must be a valid function handle>
+%!  decic (1, 0, [1; 0; 0], [1; 1; 0], [-1e-4; 1; 0], [0; 0; 0]);
+%!error <INIT must be a numeric scalar value>
+%!  decic (@rob, [1, 1], [1; 0; 0], [1; 1; 0], [-1e-4; 1; 0], [0; 0; 0]);
+%!error <length of y0, fixed_y0, yp0 and fixed_yp0 must be equal>
+%!  decic (@rob, 0, [0; 0], [1; 1; 0], [-1e-4; 1; 0], [0; 0; 0]);
+%!error <y0, fixed_y0, yp0 and fixed_yp0 must be numeric vectors>
+%!  decic (@rob, 0, [1; 0; 0], [1; 0],"", [0; 0; 0]);
+%!error <y0, fixed_y0, yp0 and fixed_yp0 must be numeric vectors>
+%!  decic (@rob, 0, [1; 0; 0], [1; 1; 0], [-1e-4; 1; 0], [0; 0; "1"]);
+%!error <fixed_y0 and fixed_yp0 must be boolean vectors>
+%!  decic (@rob, 0, [1; 0; 0], [5; 5; 0], [-1e-4; 1; 0], [0; 0; 0]);
+%!error <fixed_y0 and fixed_yp0 must be boolean vectors>
+%!  decic (@rob, 0, [1; 0; 0], [1; 1; 0], [-1e-4; 1; 0], [0; 4; 0]);
+%!error #too many components fixed
+%!  decic (@rob, 0, [1; 0; 0], [1; 1; 0], [-1e-4; 1; 0], [0; 1; 1]);
+
+

--- a/scripts/ode/module.mk	Tue Aug 23 02:13:28 2016 +0200
+++ b/scripts/ode/module.mk	Tue Aug 23 02:31:51 2016 +0200
@@ -15,11 +15,14 @@
   scripts/ode/private/starting_stepsize.m
 
 scripts_ode_FCN_FILES =  \
+  scripts/ode/ode15i.m \
+  scripts/ode/ode15s.m \
   scripts/ode/ode23.m \
   scripts/ode/ode45.m \
   scripts/ode/odeset.m \
   scripts/ode/odeget.m \
-  scripts/ode/odeplot.m
+  scripts/ode/odeplot.m \
+  scripts/ode/decic.m
 
 scripts_odedir = $(fcnfiledir)/ode
 

--- a/scripts/ode/ode23.m	Tue Aug 23 02:13:28 2016 +0200
+++ b/scripts/ode/ode23.m	Tue Aug 23 02:31:51 2016 +0200
@@ -151,6 +151,7 @@
            "ode23: FUN must be a valid function handle");
   endif
 
+
   ## Start preprocessing, have a look which options are set in odeopts,
   ## check if an invalid or unused option is set.
   [defaults, classes, attributes] = odedefaults (numel (init),