function varargout = Enzyme( x, flag, uds)
%
%  Enzyme Reaction Problem due to Kowalik/Osborne
%     
%  Used as example 4.1 in
%  P. Deuflhard: Newton Methods for Nonlinear Problems.
%  Affine Invariance and Adaptive Algorithms.
%  Series Computational Mathematics 35 
%  Springer (2004)
%
%  IN:
%    x    - the column vector of length n of parameters
%    flag - lower case string indicating the type of information 
%           that the function should return
%    uds  - a user defined structure
%
%  OUT:
%  According to the value of flag, the function must return
%    flag         varargout
%    'init'     - [fobs, fwght, x0, options, uds] 
%                 column vector of lentgth m with observations
%                 column vector of lentgth m with weights
%                 column vector of lentgth n with an initial guess for
%                 the model parameters to be estimated
%                 structure with input options for NLSQ solver, may be []
%                 user defined structure
%    ''         - [f]
%                 column vector of lentgth m with the values of the
%                 model function fmodel(x)
%    'jacobian' - [J]
%                 Jacobian (m,n)-matrix J(x) = df/dx
%
%  NOTE
%    This version uses a relative residual, i.e., fwght = fobs
%

  if nargin < 2
    flag = 'init';
  end
  if nargin < 3 | isstruct(uds) == 0
    uds = [];
  end
 
  switch flag
  case ''                           % Return model function f(x).
    varargout{1} = f(x,uds) ;
  case 'jacobian'                   % Return Jacobian matrix df/dx.
   varargout{1} = jacobian(x,uds);
  case 'init'                       % Return [fo, fw, x0, options, udf].
    [varargout{1:5}] = init(uds) ;
  otherwise
    error([' Error in function Enzyme, unknown flag ''' flag '''.']) ;
  end

 
function fx = f( x, uds)
  t = uds.ti;
  fx = (x(1).*(t.^2 + x(2)*t)) ./ (t.^2 + x(3)*t + x(4));
   
 
function J = jacobian( x, uds)
  t = uds.ti;
  u = x(1).*(t.^2 + x(2)*t);
  v = t.^2 + x(3)*t + x(4);
  J(:,1) = (t.^2 + x(2)*t) ./ v;
  J(:,2) = (x(1)*t) ./ v;
  J(:,4) = - u ./ v.^2;
  J(:,3) = J(:,4) .* t;
 

function [fobs, fwght, x0, options, uds] = init( uds)

  t = [ 4, 2, 1, 0.5, 0.25, 0.167, 0.125, 0.1, ...
       0.0823, 0.0714, 0.0625]';

  fobs = [ 0.1957, 0.1947, 0.1735, 0.1600, ...
            0.0844, 0.0627, 0.0456, 0.0342, ...
            0.0323, 0.0235, 0.0246]';
            
  m = length( fobs);   
  fwght = ones( m, 1);    % to get an absolute residual
  fwght = fobs;           % to get an relative residual
   
  x0 = [0.25, 0.39, 0.415 0.39]' ;
%  x0 = [0.19, 0.20, 0.12 0.14]' ;

  options = [] ;
  options.xthrsh = 0.1*ones( 4, 1);  
  options.numdif = 0;
%
  uds.ti = t;
%     End of problem definition for Enzyme Reaction Problem


