function varargout = problem(x,flag,p)
%
%     Chebyquad function.
%
%     VAR = PROBLEM(X,FLAG,PAR) computes the problem function f(X),
%       Jacobian df/dx(X) and startvector for the Newton iteration
%       of the problem.
%
%     FX = PROBLEM(X,'',par)          
%          returns the right hand side f(X) for a input column vector X.
%     JAC = PROBLEM(X,'jacobian',par) 
%           returns the Jacobian df/dx(X) for a input column vector X.
%     [X0,OPTIONS] = PROBLEM(X,'init',par)
%        returns a startvector X0 and problem specific solver options
%        OPTIONS which should override any solver default options.
%        For more help on OPTIONS, see help nleq1.
%
%     Valid function parameters PAR are:
%     PAR.n : the dimension n of the Chebyquad problem -
%             n may be range from 2 to 9. 
%               Default is 2.
%             Note that for n=8 there exists no solution of the problem.
%
 switch flag
 case ''                                 % Return y = f(x).
   [varargout{1:2}] = f(x) ;
 case 'jacobian'                         % Return Jacobian matrix df/dx.
   [varargout{1:2}] = jacobian(x);
 case 'init'                             % Return default [x0,options].
   [varargout{1:2}] = init(x,p) ;
 otherwise
   error(['Unknown flag ''' flag '''.']) ;
 end

% --------------------------------------------------------------------------
 
function [fx,fail] = f(x)
   fail = 0;
   [n,dummy] = size(x) ;
   fx = zeros(n,1) ;
   for i=2:2:n
	 fx(i)=n/(i*i-1) ;
   end
   for l=1:n
	 factt = 4.0*x(l)-2.0 ;
	 ti2 = 1.0 ;
	 ti1 = 0.5*factt ;
	 fx(1)=ti1+fx(1) ;
	 for i=2:n
	   ti = factt*ti1-ti2 ;
	   fx(i)=ti+fx(i) ;
	   ti2 = ti1 ;
	   ti1 = ti ;
	 end
   end

% --------------------------------------------------------------------------
 
function [J,fail] = jacobian(x)
   fail = 0;
   [n,dummy] = size(x) ;
   J = zeros(n) ;
   for j=1:n
	 factt = 4.0*x(j)-2.0 ;
	 ti2 = 1.0 ;
	 ti1 = 0.5*factt ;
	 tabli2 = 0.0 ;
	 tabli1 = 2.0 ;
	 J(1,j)=tabli1 ;
	 for i=2:n
	   ti = factt*ti1-ti2 ;
	   tabli = 4.0*ti1+factt*tabli1-tabli2 ;
	   J(i,j)=tabli ;
	   ti2 = ti1 ;
	   ti1 = ti ;
	   tabli2 = tabli1 ;
	   tabli1 = tabli ;
	 end
   end


% --------------------------------------------------------------------------
 
function [x0,options] = init(x,p)
  if isempty(p)
    n = 2 ;
  else
    n = p.n ;
  end
  x0 = [1:n]' / (n+1) ;
  options = [] ;
  
%     End of problem definition Chebyquad problem
  

