function C=SolveEquation(TM,cg,method)
% SolveEquation - solve tensor-matrix equation
%
% Syntax: C=SolveEquation(TM,cg,method)
%         C=SolveEquation(TM,cg)
%         C=SolveEquation(TM)
%
% TM must be a TensorMatrix object (a discretized PDE).
% Peclet is the Peclet number (surprise).
%
% C returned is a matrix (or, if cg is supplied, a NumField over cg)
% such that 
%
%     A_ijkl C_kl = B_ij
%
% (Here we use the Einstein convention for summation over 
% repeated indeces)
%
% Bugs and Features: Currently this is done with Gaussian 
% elimination. See 'help sparfun' for alternatives which 
% could be used in stead.
%
% See also: TensorMatrix, NumField, compgrid
%

if nargin<3,
	C = TM.A\TM.B;
else
	C = mldivide(TM.A,TM.B,method);
end

if (nargin>1) & (~isempty(cg)),
   C = NumField(C,cg);   
end