function R = solve(self,tlist,options)
% solve - Solve the thermal model
% R = solve(thermalmodel) solves a steady-state thermal model and
% returns the result in a SteadyStateThermalResults object.
%
% R = solve(thermalmodel, tlist) solves a transient thermal model and
% returns the result in a TransientThermalResults object. The input
% argument tlist is a monotonic increasing or decreasing vector that
% specifies the times at which the solution is requested.
%
% R = solve(thermalmodel,"DecayRange",[lambda1,lambda2]) solves a
% thermal modal analysis model for all modes in the decay range
% [lambda1,lambda2] and returns the result in a ModalThermalResults
% object.

% Example 1: Solve a steady-state thermal model.
%
%     %Create a thermal model for this problem.
%     thermalmodel = createpde('thermal');
%
%     %Import and plot the block geometry.
%     importGeometry(thermalmodel,'Block.stl');
%
%     %Assign material properties.
%     thermalProperties(thermalmodel,'ThermalConductivity',0.80);
%
%     % Apply a constant temperature of 100 to the left side of the block
%     %(face 1) and a constant temperature of 300 at the right side of the
%     %block (face 3). All other faces are insulated by default.
%     thermalBC(thermalmodel,'Face',1,'Temperature',100);
%     thermalBC(thermalmodel,'Face',3,'Temperature',300);
%
%     %Mesh the geometry and solve the problem.
%     generateMesh(thermalmodel);
%     thermalresults = solve(thermalmodel)
%
%     %Plot temperature at nodal locations.
%     pdeplot3D(thermalmodel,'ColorMapData',thermalresults.Temperature)
%
%  Example 2: Solve transient thermal model.
%
%     thermalmodel = createpde('thermal','transient');
%     gm = @squareg;
%     geometryFromEdges(thermalmodel,gm);
%     generateMesh(thermalmodel);
%     thermalProperties(thermalmodel,'ThermalConductivity',50,...
%                                     'MassDensity',2500,...
%                                     'SpecificHeat',600);
%     bc = thermalBC(thermalmodel,'Edge',[1,3],'Temperature',100);
%     qs = internalHeatSource(thermalmodel,2);
%     ic = thermalIC(thermalmodel,0);
%     tlist = 0:0.1:1;
%     thermalresults = solve(thermalmodel,tlist)
%     pdeplot(thermalmodel,'xydata',thermalresults.Temperature)
%
% See also pde.ThermalModel, pde.ThermalModel/thermalBC,
%           pde.ThermalModel/thermalIC, pde.ThermalModel/internalHeatSource
%           pde.ThermalModel/thermalProperties

% Copyright 2016-2023 The MathWorks, Inc.

arguments
    self pde.ThermalModel
    tlist (1,:) = [] % Optional positional argument.
    options.DecayRange (1,:) {mustBeA(options.DecayRange,"double"),mustBeReal,mustBeVector,mustBeNonsparse,mustBeNonNan,mustBeNonempty};
    options.ModalResults pde.ModalThermalResults;
    options.Snapshots {mustBeA(options.Snapshots,"double"),mustBeReal,mustBeNonsparse,mustBeNonNan,mustBeFinite,mustBeNonempty}; 
end

narginchk(1,4);
checkics = false;
if self.IsTimeDependent
    checkics = true;
    validateattributes(tlist,{'numeric'},{'real', 'finite', 'nonsparse', 'nonnan'});
end
performSolverPrecheck(self, checkics);

if nargin == 2 && ~self.IsTimeDependent
    warning(message('pde:pdeModel:stationaryTlistIgnored'))
elseif nargin == 1 && self.IsTimeDependent
    error(message('pde:pdeModel:timeDependentNoTlist'))
end

allbcsnumeric = true;
if ~isempty(self.BoundaryConditions)
    numbcs = numel(self.BoundaryConditions.ThermalBCAssignments);
    for i = 1:numbcs
        thisbc = self.BoundaryConditions.ThermalBCAssignments(i);
        if ~thisbc.numericCoefficients() || ~isempty(thisbc.Emissivity)
            allbcsnumeric = false;
            break;
        end
    end
end

[coefstruct,constCoeffs] = packCoefficients(self);
constCoeffsBCs = constCoeffs && allbcsnumeric;

if contains(self.AnalysisType,"modal")
    if ~isfield(options,'DecayRange') && ~isfield(options,'Snapshots')
        error(message('pde:StructuralModel:noFreqRangeSpecified','DecayRange or Snapshots'));
    end
    
    Tmean = [];
    if isfield(options,'DecayRange')
        tsecondorder = false;
        [Phi,lambda,B] = solveEigenvalue(self,options.DecayRange,coefstruct,tsecondorder);
        modeType = "EigenModes";
    else
        [Phi,lambda,B,Tmean] = solvePODModes(self,options.Snapshots,coefstruct);
        modeType = "PODModes";
    end
   
    if isempty(lambda)
        R = pde.ModalThermalResults;
    else
        R = pde.ModalThermalResults(self,Phi,lambda,B,modeType,Tmean);
    end
    return
end

if self.IsTimeDependent
    self.checkTlist(tlist);
    u0 = self.InitialConditions.packInitialConditions();
    if ~isfield(options,'ModalResults')
        % Direct solve
        u = self.solveTimeDependent(coefstruct,u0,[],tlist,false);
    else
        % Modal solver
        u = self.solveThermalModalTimeDependent(coefstruct,u0,tlist,options.ModalResults);
    end
    R = pde.TransientThermalResults(self,u,tlist);
else
    [p,e,t] = self.Mesh.meshToPet();
    if ~isempty(self.InitialConditions)
        u0 = self.InitialConditions.packInitialConditions();
    else
        u0 = pdeuxpd(p,0,self.PDESystemSize);
    end
    femodel = pde.DiscretizedPDEModel(self,p,e,t,coefstruct,u0,constCoeffsBCs);
    if (~femodel.IsStationaryNonlinear)
        if ~isempty(self.InitialConditions)
            warning(message('pde:ThermalModel:noICforSteadyStateLinearModel'));
        end
        u = self.solveStationary(femodel);
        R = pde.SteadyStateThermalResults(self,u);
    else
        u = self.solveStationaryNonlinear(coefstruct, u0);
        R = pde.SteadyStateThermalResults(self,u);
    end
end


end

