function [ x, fval, output ] = fminode ( x0, m );
% FMinODE
%
% This m-file is part of the material of the lectures on
%   Numerical Methods in Nonlinear Optimization
% at the Math Department of Halle University, Germany.
%
% Author :      Martin Arnold, martin.arnold@mathematik.uni-halle.de
% Version of :  Apr 25, 2006
%
% Parameter estimation involving ODE initial value problems
%
% Parameters:
%   x0     (input)  : initial guess  x0 = [ a, b ]
%   m      (input)  : number of measured data
%   x      (output) : optimal value computed by fminsearch
%   fval   (output) : optimal function value
%   output (output) : structure with statistical data of fminsearch
%
% Example:
%   [ x, fval, output ] = fminode ( [  0.4, 2.0 ], 201 )
%   [ x, fval, output ] = fminode ( [  1.0, 5.0 ], 101 )
%   [ x, fval, output ] = fminode ( [ -1.0, 3.0 ], 401 )

% -> generate test data
a  = 0.4;
b  = 2.0;
te = 1.0;
m  = 201;

x = [ a, b ];

y0 = b;
ts = linspace ( 0, te, m );

rhs = @(t,y)( tdglrhs(t,y,x) );
options = odeset ( 'RelTol', 1.0e-8, 'AbsTol', 1.0e-8 );

[ t, y ] = ode15s ( rhs, ts, y0, options );

eta = y(2:end);

% -> problem and solver setup
mytdglf = @(x) ( tdglf ( x, m, ts, eta ) );

options = optimset (                            ...
                     'Display',     'off'       ...
                   , 'TolX',        1e-6        ...
                   , 'MaxFunEvals', 5000        ...
                   , 'MaxIter',     2000        ...
                   );

% -> minimization
[ x, fval, exitflag, output ] = fminsearch ( mytdglf, x0, options );

% ----------------------------------------------------------------- %

function f = tdglf ( x, m, ts, eta );
% tDglF
%
% Function to be minimized
%
% Paramters:
%   x      (input)  : current iterate
%   m      (input)  : number of data points
%   ts     (input)  : points of measurements
%   eta    (input)  : measured data
%   f      (output) : function value

% -> ODE time integration
y0 = x(end);

rhs = @(t,y)( tdglrhs(t,y,x) );
options = odeset ( 'RelTol', 1.0e-8, 'AbsTol', 1.0e-8 );

[ t, y ] = ode15s ( rhs, ts, y0, options );

% -> evaluate f(x)
f = sum ( (y(2:end)-eta).^2 );

% ----------------------------------------------------------------- %

function phi = tdglrhs ( t, y, x );
% tDglRhs
%
% ODE right hand side
%
% Paramters:
%   t      (input)  : current time
%   y      (input)  : current state
%   phi    (output) : current derivative

% -> read parameter
a = x(1);

% -> evaluate right hand side
phi = a * y;