[fnc,neq,x,rtolx,atolx,tolf,iflag,rw,lrw,iw,liw]=sos(fnc,neq,x,rtolx,atolx,tolf,iflag,rw,lrw,iw,liw);

function [fnc,neq,x,rtolx,atolx,tolf,iflag,rw,lrw,iw,liw]=sos(fnc,neq,x,rtolx,atolx,tolf,iflag,rw,lrw,iw,liw);
persistent inpflg iprint k1 k2 k3 k4 k5 k6 mxit nc ncjs nsri nsrrc xern1 xern3 xern4 ; 

if isempty(inpflg), inpflg=0; end;
if isempty(iprint), iprint=0; end;
if isempty(k1), k1=0; end;
if isempty(k2), k2=0; end;
if isempty(k3), k3=0; end;
if isempty(k4), k4=0; end;
if isempty(k5), k5=0; end;
if isempty(k6), k6=0; end;
if isempty(mxit), mxit=0; end;
if isempty(nc), nc=0; end;
if isempty(ncjs), ncjs=0; end;
if isempty(nsri), nsri=0; end;
if isempty(nsrrc), nsrrc=0; end;
%***BEGIN PROLOGUE  SOS
%***PURPOSE  Solve a square system of nonlinear equations.
%***LIBRARY   SLATEC
%***CATEGORY  F2A
%***TYPE      SINGLE PRECISION (SOS-S, DSOS-D)
%***KEYWORDS  BROWN'S METHOD, NEWTON'S METHOD, NONLINEAR EQUATIONS,
%             ROOTS, SOLUTIONS
%***AUTHOR  Watts, H. A., (SNLA)
%***DESCRIPTION
%
%     SOS solves a system of NEQ simultaneous nonlinear equations in
%     NEQ unknowns.  That is, it solves the problem   F(X)=0
%     where X is a vector with components  X(1),...,X(NEQ)  and  F
%     is a vector of nonlinear functions.  Each equation is of the form
%
%               F (X(1),...,X(NEQ))=0     for K=1,...,NEQ.
%                K
%
%     The algorithm is based on an iterative method which is a
%     variation of Newton's method using Gaussian elimination
%     in a manner similar to the Gauss-Seidel process.  Convergence
%     is roughly quadratic.  All partial derivatives required by
%     the algorithm are approximated by first difference quotients.
%     The convergence behavior of this code is affected by the
%     ordering of the equations, and it is advantageous to place linear
%     and mildly nonlinear equations first in the ordering.
%
%     Actually, SOS is merely an interfacing routine for
%     calling subroutine SOSEQS which embodies the solution
%     algorithm.  The purpose of this is to add greater
%     flexibility and ease of use for the prospective user.
%
%     SOSEQS calls the accompanying routine SOSSOL, which solves special
%     triangular linear systems by back-substitution.
%
%     The user must supply a function subprogram which evaluates the
%     K-th equation only (K specified by SOSEQS) for each call
%     to the subprogram.
%
%     SOS represents an implementation of the mathematical algorithm
%     described in the references below.  It is a modification of the
%     code SOSNLE written by H. A. Watts in 1973.
%
% **********************************************************************
%   -Input-
%
%     FNC -Name of the function program which evaluates the equations.
%          This name must be in an EXTERNAL statement in the calling
%          program.  The user must supply FNC in the form FNC(X,K),
%          where X is the solution vector (which must be dimensioned
%          in FNC) and FNC returns the value of the K-th function.
%
%     NEQ -Number of equations to be solved.
%
%     X   -Solution vector.  Initial guesses must be supplied.
%
%     RTOLX -Relative error tolerance used in the convergence criteria.
%          Each solution component X(I) is checked by an accuracy test
%          of the form   ABS(X(I)-XOLD(I)) .LE. RTOLX*ABS(X(I))+ATOLX,
%          where XOLD(I) represents the previous iteration value.
%          RTOLX must be non-negative.
%
%     ATOLX -Absolute error tolerance used in the convergence criteria.
%          ATOLX must be non-negative.  If the user suspects some
%          solution component may be zero, he should set ATOLX to an
%          appropriate (depends on the scale of the remaining variables)
%          positive value for better efficiency.
%
%     TOLF -Residual error tolerance used in the convergence criteria.
%          Convergence will be indicated if all residuals (values of the
%          functions or equations) are not bigger than TOLF in
%          magnitude.  Note that extreme care must be given in assigning
%          an appropriate value for TOLF because this convergence test
%          is dependent on the scaling of the equations.  An
%          inappropriate value can cause premature termination of the
%          iteration process.
%
%     IFLAG -Optional input indicator.  You must set  IFLAG=-1  if you
%          want to use any of the optional input items listed below.
%          Otherwise set it to zero.
%
%     RW  -A REAL work array which is split apart by SOS and used
%          internally by SOSEQS.
%
%     LRW -Dimension of the RW array.  LRW must be at least
%                    1 + 6*NEQ + NEQ*(NEQ+1)/2
%
%     IW  -An INTEGER work array which is split apart by SOS and used
%          internally by SOSEQS.
%
%     LIW -Dimension of the IW array. LIW must be at least  3 + NEQ.
%
%   -Optional Input-
%
%     IW(1) -Internal printing parameter.  You must set  IW(1)=-1  if
%          you want the intermediate solution iterates to be printed.
%
%     IW(2) -Iteration limit.  The maximum number of allowable
%          iterations can be specified, if desired.  To override the
%          default value of 50, set IW(2) to the number wanted.
%
%     Remember, if you tell the code that you are using one of the
%               options (by setting IFLAG=-1), you must supply values
%               for both IW(1) and IW(2).
%
% **********************************************************************
%   -Output-
%
%     X   -Solution vector.
%
%     IFLAG -Status indicator
%
%                         *** Convergence to a Solution ***
%
%          1 Means satisfactory convergence to a solution was achieved.
%            Each solution component X(I) satisfies the error tolerance
%            test   ABS(X(I)-XOLD(I)) .LE. RTOLX*ABS(X(I))+ATOLX.
%
%          2 Means procedure converged to a solution such that all
%            residuals are at most TOLF in magnitude,
%            ABS(FNC(X,I)) .LE. TOLF.
%
%          3 Means that conditions for both IFLAG=1 and IFLAG=2 hold.
%
%          4 Means possible numerical convergence.  Behavior indicates
%            limiting precision calculations as a result of user asking
%            for too much accuracy or else convergence is very slow.
%            Residual norms and solution increment norms have
%            remained roughly constant over several consecutive
%            iterations.
%
%                         *** Task Interrupted ***
%
%          5 Means the allowable number of iterations has been met
%            without obtaining a solution to the specified accuracy.
%            Very slow convergence may be indicated.  Examine the
%            approximate solution returned and see if the error
%            tolerances seem appropriate.
%
%          6 Means the allowable number of iterations has been met and
%            the iterative process does not appear to be converging.
%            A local minimum may have been encountered or there may be
%            limiting precision difficulties.
%
%          7 Means that the iterative scheme appears to be diverging.
%            Residual norms and solution increment norms have
%            increased over several consecutive iterations.
%
%                         *** Task Cannot Be Continued ***
%
%          8 Means that a Jacobian-related matrix was singular.
%
%          9 Means improper input parameters.
%
%          *** IFLAG should be examined after each call to   ***
%          *** SOS with the appropriate action being taken.  ***
%
%
%     RW(1) -Contains a norm of the residual.
%
%     IW(3) -Contains the number of iterations used by the process.
%
% **********************************************************************
%***REFERENCES  K. M. Brown, Solution of simultaneous nonlinear
%                 equations, Algorithm 316, Communications of the
%                 A.C.M. 10, (1967), pp. 728-729.
%               K. M. Brown, A quadratically convergent Newton-like
%                 method based upon Gaussian elimination, SIAM Journal
%                 on Numerical Analysis 6, (1969), pp. 560-569.
%***ROUTINES CALLED  SOSEQS, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   801001  DATE WRITTEN
%   890831  Modified array declarations.  (WRB)
%   890831  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900510  Convert XERRWV calls to XERMSG calls, changed Prologue
%           comments to agree with DSOS.  (RWC)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  SOS
x_shape=size(x);x=reshape(x,1,[]);
rw_shape=size(rw);rw=reshape(rw,1,[]);
iw_shape=size(iw);iw=reshape(iw,1,[]);
if isempty(xern1), xern1=repmat(' ',1,8); end;
if isempty(xern3), xern3=repmat(' ',1,16); end;
if isempty(xern4), xern4=repmat(' ',1,16); end;
%***FIRST EXECUTABLE STATEMENT  SOS
inpflg = fix(iflag);
%
%     CHECK FOR VALID INPUT
%
if( neq<=0 )
xern1=sprintf(['%8i'], neq);
xermsg('SLATEC','SOS',['THE NUMBER OF EQUATIONS ',['MUST BE A POSITIVE INTEGER.  YOU HAVE CALLED THE ',['CODE WITH NEQ = ',xern1]]],1,1);
iflag = 9;
end;
%
if( rtolx<0.0d0 || atolx<0.0d0 )
xern3=sprintf([repmat('%15.6f',1,1)], atolx);
xern4=sprintf([repmat('%15.6f',1,1)], rtolx);
xermsg('SLATEC','SOS',['THE ERROR TOLERANCES FOR ',['THE SOLUTION ITERATES CANNOT BE NEGATIVE. YOU HAVE ',['CALLED THE CODE WITH  RTOLX = ',[xern3,[' AND ATOLX = ',xern4]]]]],2,1);
iflag = 9;
end;
%
if( tolf<0.0d0 )
xern3=sprintf([repmat('%15.6f',1,1)], tolf);
xermsg('SLATEC','SOS',['THE RESIDUAL ERROR ',['TOLERANCE MUST BE NON-NEGATIVE.  YOU HAVE CALLED THE ',['CODE WITH TOLF = ',xern3]]],3,1);
iflag = 9;
end;
%
iprint = 0;
mxit = 50;
if( inpflg==(-1) )
if( iw(1)==(-1) )
iprint = -1;
end;
mxit = fix(iw(2));
if( mxit<=0 )
xern1=sprintf(['%8i'], mxit);
xermsg('SLATEC','SOS',['YOU HAVE TOLD THE CODE ',['TO USE OPTIONAL IN PUT ITEMS BY SETTING  IFLAG=-1. ',['HOWEVER YOU HAVE CALLED THE CODE WITH THE MAXIMUM ',['ALLOWABLE NUMBER OF ITERATIONS SET TO  IW(2) = ',xern1]]]],4,1);
iflag = 9;
end;
end;
%
nc =fix(fix((neq.*(neq+1))./2));
if( lrw<1+6.*neq+nc )
xern1=sprintf(['%8i'], lrw);
xermsg('SLATEC','SOS',['DIMENSION OF THE RW ARRAY ',['MUST BE AT LEAST 1 + 6*NEQ + NEQ*(NEQ+1)/2 .  YOU HAVE ',['CALLED THE CODE WITH LRW = ',xern1]]],5,1);
iflag = 9;
end;
%
if( liw<3+neq )
xern1=sprintf(['%8i'], liw);
xermsg('SLATEC','SOS',['DIMENSION OF THE IW ARRAY ',['MUST BE AT LEAST   3 + NEQ.  YOU HAVE CALLED THE CODE ',['WITH  LIW = ',xern1]]],6,1);
iflag = 9;
end;
%
if( iflag~=9 )
ncjs = 6;
nsrrc = 4;
nsri = 5;
%
k1 = fix(nc + 2);
k2 = fix(k1 + neq);
k3 = fix(k2 + neq);
k4 = fix(k3 + neq);
k5 = fix(k4 + neq);
k6 = fix(k5 + neq);
%
[fnc,neq,x,rtolx,atolx,tolf,iflag,mxit,ncjs,nsrrc,nsri,iprint,dumvar13,dumvar14,nc,dumvar16,dumvar17,dumvar18,dumvar19,dumvar20,dumvar21,iw(sub2ind(size(iw),max(4,1)):end)]=soseqs(fnc,neq,x,rtolx,atolx,tolf,iflag,mxit,ncjs,nsrrc,nsri,iprint,rw(1),rw(2:2+nc-1),nc,rw(sub2ind(size(rw),max(k1,1)):end),rw(sub2ind(size(rw),max(k2,1)):end),rw(sub2ind(size(rw),max(k3,1)):end),rw(sub2ind(size(rw),max(k4,1)):end),rw(sub2ind(size(rw),max(k5,1)):end),rw(sub2ind(size(rw),max(k6,1)):end),iw(sub2ind(size(iw),max(4,1)):end));   dumvar13i=find((rw(1))~=(dumvar13));dumvar14i=find((rw(2:2+nc-1))~=(dumvar14));dumvar16i=find((rw(sub2ind(size(rw),max(k1,1)):end))~=(dumvar16));dumvar17i=find((rw(sub2ind(size(rw),max(k2,1)):end))~=(dumvar17));dumvar18i=find((rw(sub2ind(size(rw),max(k3,1)):end))~=(dumvar18));dumvar19i=find((rw(sub2ind(size(rw),max(k4,1)):end))~=(dumvar19));dumvar20i=find((rw(sub2ind(size(rw),max(k5,1)):end))~=(dumvar20));dumvar21i=find((rw(sub2ind(size(rw),max(k6,1)):end))~=(dumvar21));   rw(1-1+dumvar13i)=dumvar13(dumvar13i); rw(2-1+dumvar14i)=dumvar14(dumvar14i); rw(k1-1+dumvar16i)=dumvar16(dumvar16i); rw(k2-1+dumvar17i)=dumvar17(dumvar17i); rw(k3-1+dumvar18i)=dumvar18(dumvar18i); rw(k4-1+dumvar19i)=dumvar19(dumvar19i); rw(k5-1+dumvar20i)=dumvar20(dumvar20i); rw(k6-1+dumvar21i)=dumvar21(dumvar21i); 
%
iw(3) = fix(mxit);
end;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
rw_shape=zeros(rw_shape);rw_shape(:)=rw(1:numel(rw_shape));rw=rw_shape;
iw_shape=zeros(iw_shape);iw_shape(:)=iw(1:numel(iw_shape));iw=iw_shape;
end
%DECK SOSSOL