function varargout = solve(varargin)
%求各種類型方程組的解析解, 
%當找不到解析解時,將自動尋求原點附近的一個近似解, 并可達任意精度.
%用法
%   solve('方程1','方程2',...,'方程N','變量1','變量2',...,'變量M')
%SOLVE  Symbolic solution of algebraic equations.
%   SOLVE('eqn1','eqn2',...,'eqnN')
%   SOLVE('eqn1','eqn2',...,'eqnN','var1,var2,...,varN')
%   SOLVE('eqn1','eqn2',...,'eqnN','var1','var2',...'varN')
%
%   The eqns are symbolic expressions or strings specifying equations.  The
%   vars are symbolic variables or strings specifying the unknown variables.
%   SOLVE seeks zeros of the expressions or solutions of the equations.
%   If not specified, the unknowns in the system are determined by FINDSYM.
%   If no analytical solution is found and the number of equations equals
%   the number of dependent variables, a numeric solution is attempted.
%
%   Three different types of output are possible.  For one equation and one
%   output, the resulting solution is returned, with multiple solutions to
%   a nonlinear equation in a symbolic vector.  For several equations and
%   an equal number of outputs, the results are sorted in lexicographic
%   order and assigned to the outputs.  For several equations and a single
%   output, a structure containing the solutions is returned.
%
%   Examples:
%
%      solve('p*sin(x) = r') chooses 'x' as the unknown and returns
%
%        ans =
%        asin(r/p)
%
%      [x,y] = solve('x^2 + x*y + y = 3','x^2 - 4*x + 3 = 0') returns
%
%        x =
%        [ 1]
%        [ 3]
%
%        y =
%        [    1]
%        [ -3/2]
%
%      S = solve('x^2*y^2 - 2*x - 1 = 0','x^2 - y^2 - 1 = 0') returns
%      the solutions in a structure.
%
%        S =
%          x: [8x1 sym]
%          y: [8x1 sym]
%
%      [u,v] = solve('a*u^2 + v^2 = 0','u - v = 1') regards 'a' as a
%      parameter and solves the two equations for u and v.
%
%      S = solve('a*u^2 + v^2','u - v = 1','a,u') regards 'v' as a
%      parameter, solves the two equations, and returns S.a and S.u.
%
%      [a,u,v] = solve('a*u^2 + v^2','u - v = 1','a^2 - 5*a + 6') solves
%      the three equations for a, u and v.
%
%      [x,y] = solve('sin(x+y)-exp(x)*y = 0','x^2-y = 2') cannot find
%      an analytic solution, so returns a numeric solution.
%
%   See also DSOLVE.

%   Copyright (c) 1993-98 by The MathWorks, Inc.
%   $Revision: 1.13 $  $Date: 1997/11/29 01:06:35 $

% Set the Explicit solution (for equations of order 4 or less)
% option on in the Maple Workspace.

maple('_EnvExplicit := true;');

% Collect input arguments together in either equation or variable lists.

eqns = [];
vars = [];
for k = 1:nargin
   v = varargin{k};
   if all(isletter(v) | ('0'<=v & v<='9') | v == '_' | v == ',')
      vars = [vars ',' v];
   else
      [t,stat] = maple(v);
      if stat
         error(['''' v ''' is not a valid expression or equation.'])
      end
      eqns = [eqns ',' v];   
   end
end

if isempty(eqns)
   warning('List of equations is empty.')
   varargout = cell(1,nargout);
   varargout{1} = sym([]);
   return
else
   eqns(1) = '{'; eqns(end+1) = '}';
end
neqns = sum(commas(eqns)) + 1;

% Determine default variables and sort variable list.

if isempty(vars)
   vars = ['[' findsym(sym(eqns),neqns) ']'];
else
   vars(1) = '['; vars(end+1) = ']';
end
v = vars;
[vars,stat] = maple('sort',v);
if stat
   error(['''' v ''' is not a valid variable list.'])
end
vars(1) = '{'; vars(end) = '}';
nvars = sum(commas(vars)) + 1;

if neqns ~= nvars
   warning([int2str(neqns) ' equations in ' int2str(nvars) ' variables.'])
end

% Seek analytic solution.

[R,stat] = maple('solve',eqns,vars);

% If no analytic solution, seek numerical solution.

if (isempty(R) & (nvars == neqns))
   [R,stat] = maple('fsolve',eqns,vars);
end

if stat
   error(R)
end

% If still no solution, give up.

if isempty(R) | findstr(R,'fsolve')
   warning('Explicit solution could not be found.');
   varargout = cell(1,nargout);
   varargout{1} = sym([]);
   return
end

% Expand any RootOf.

while ~isempty(findstr(R,'RootOf'))
   k = min(findstr(R,'RootOf'));
   p = findstr(R,'{'); p = max(p(p<k));
   q = findstr(R,'}'); q = min(q(q>k));
   s = R(p:q);
   t = s(min(findstr(s,'RootOf'))+6:end);
   e = min(find(cumsum((t=='(')-(t==')'))==0));
   if isempty(find(commas(t(2:e-1))))
      % RootOf with one argument, possibly an imbedded RootOf.
      [s,stat] = maple('allvalues',s,'dependent');
      if stat
         error(s)
      end
   else
      % RootOf with a good estimate as the second argument.
      s = maple('evalf',s);
   end
   R = [R(1:p-1) s R(q+1:end)];
end

% Parse the result.

if nvars == 1 & nargout <= 1

   % One variable and at most one output.
   % Return a single scalar or vector sym.

   S = sym([]);
   c = find(commas(R) | R == '}');
   for p = find(R == '=')
      q = min(c(c>p));
      t = trim(R(p+1:q-1));  % The solution (xk)
      S = [S; sym(t)];
   end
   varargout{1} = S;

else

   % Several variables.
   % Create a skeleton structure.

   c = [1 find(commas(vars)) length(vars)];
   S = [];
   for j = 1:nvars
      v = trim(vars(c(j)+1:c(j+1)-1));
      S = setfield(S,v,[]);
   end

   % Complete the structure.

   c = [1 find(commas(R) | R == '{' | R == '}') length(R)];
   for p = find(R == '=')
      q = max(c(c<p));
      v = trim(R(q+1:p-1));  % The variable (x)
      q = min(c(c>p));
      t = trim(R(p+1:q-1));  % The solution (xk)
      S = setfield(S,v,[getfield(S,v); sym(t)]);
   end
   
   if nargout <= 1

      % At most one output, return the structure.
      varargout{1} = S;

   elseif nargout == nvars

      % Same number of outputs as variables.
      % Match results in lexicographic order to outputs.
      v = fieldnames(S);
      for j = 1:nvars
         varargout{j} = getfield(S,v{j});
      end

   else
      error([int2str(nvars) ' variables does not match ' ...
             int2str(nargout) ' outputs.'])
   end
end

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

function s = trim(s);
%TRIM  TRIM(s) deletes any leading or trailing blanks.
while s(1) == ' ', s(1) = []; end
while s(end) == ' ', s(end) = []; end

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

function c = commas(s)
%COMMAS  COMMAS(s) is true for commas not inside parentheses.
p = cumsum((s == '(') - (s == ')'));
c = (s == ',') & (p == 0);