function sol = integrateObjDisc(m, con, obj, opts)
% m is a scalar
% con is a scalar
% obj is a vector

% Constants
nX = m.nX;
nObj = numel(obj);

% Construct system
[gDer, gJac, gEvents] = constructSystem();

%TODO: allow for starting from steady-state

% Initial conditions
if opts.UseModelICs
    ic = [m.ic; 0];
else
    ic = [con.ic; 0];
end

% Get the times that the solver will need to stop
stopTimes = unique([0; con.tF; con.tStop((con.tStop < con.tF))]);

% Integrate
sol = accumulateSol(gDer, gJac, gEvents, [], ic, con, opts.AbsTol(1:nX+1), opts.RelTol, stopTimes, 1, 1:m.nX, opts.tGet);

% End of function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% The system for integrating x and g %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [gDer, gJac, gEvents] = constructSystem()
        
        f      = m.f;
        dfdx   = m.dfdx;
        
        gDer = @der;
        gJac = @jac;
        if ~isempty(con.events)
            events  = con.events.events;
            gEvents = @eventFun;
        else
            gEvents = [];
        end
        
        % Derivative of [x; G] with respect to time
        function val = der(t, joint, u)
            u = u(t);            
            x = joint(1:nX);
            
            % Sum continuous objective functions
            g = 0;
            for i = 1:nObj
                g = g + opts.ObjWeights(i) * obj(i).g(t, x, u);
            end
            
            val = [f(t, x, u); g];
        end
        
        % Jacobian of [x; G] derivative
        function val = jac(t, joint, u)
            u = u(t);
            x = joint(1:nX);
            
            % Sum continuous objective gradients
            dgdx = zeros(1,nX);
            for i = 1:nObj
                dgdx = dgdx + opts.ObjWeights(i) * obj(i).dgdx(t, x, u);
            end
            
            val = [dfdx(t, x, u), sparse(nX,1);
                            dgdx,            0];
        end
        
        % Event function for system
        function [val, dir, term] = eventFun(t, joint, u)
            u = u(t);
            x = joint(1:nX);
            [val, dir, term] = events(t, x, u);
        end
    end
end