% Copyright (C) 2001 Michel Juillard
%
function [x,check] = dynare_solve(func,x,varargin)
  global gstep_
  
  check = 0;
  func = str2func(func);
  if ~isempty(which('fsolve')) & sscanf(version('-release'),'%d') >= 13;
    options=optimset('fsolve');
    options.MaxFunEvals = 20000;
    options.Display = 'off';
    [x,fval,exitval,output] = fsolve(func,x,options,varargin{:});
    if exitval > 0
      check = 0;
    else
      check = 1;
    end
    return
  end

  nn = size(x,1) ;
  tolf = eps^(2/3) ;

  fjac = zeros(nn,nn) ;

  fvec = feval(func,x,varargin{:});

  i = find(~isfinite(fvec));
  
  if ~isempty(i)
    disp(['STEADY:  numerical initial values incompatible with the following' ...
	  ' equations'])
    disp(i')
    error('exiting ...')
  end
  
  f = 0.5*fvec'*fvec ;

  if max(abs(fvec)) < 0.01*tolf
    return ;
  end

  dh = max(abs(x),gstep_*ones(nn,1))*eps^(1/3);
  for j = 1:nn
    xdh = x ;
    xdh(j) = xdh(j)+dh(j) ;
    fjac(:,j) = (feval(func,xdh,varargin{:}) - fvec)./dh(j) ;
  end

  [j1,j2,r,s] = dmperm(fjac);
  
  for i=length(r)-1:-1:1
    [x,check]=solve1(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),varargin{:});
    if check
      error(sprintf('Solve block = %d check = %d\n',i,check));
    end
  end
  
