function sol = integrateObjSens(m, con, obj, opts)

% Constants
nX = m.nX;
nVP = sum(sum(opts.UseParams));
nVX = sum(sum(opts.UseICs));
nV  = nVP + nVX;
nObj = numel(obj);

% Construct system
[jDer, jJac, jEvents] = constructSystem();

%TODO: allow for starting from steady-state

% Initial conditions [x0; G0; vec(dxdv0); vec(dGdv0)]
if opts.UseModelICs
    x0 = m.ic;
else
    x0 = con.ic;
end

% Initial effect of rates on sensitivities is 0
dxdp0 = zeros(nX, nVP); % Active rate parameters

% Initial effect of ics on sensitivities is 1 for that species
dxdx0                = zeros(nX,nVX);
dxdx0(opts.UseICs,:) = eye(nVX);

% Combine them into a vector
j0 = [x0; 0; vec([dxdp0, dxdx0]); zeros(nV,1)];

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

% Integrate [x; G; dxdv; dGdv] with respect to time
sol = accumulateSol(jDer, jJac, jEvents, [], j0, con, opts.AbsTol(1:nX+1+nX*nV+nV), opts.RelTol, stopTimes, 1, 1:nX);

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

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% The system for integrating x, G, dxdp, and D %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [jDer, jJac, jEvents] = constructSystem()
        
        Ip = speye(nV);

        dxdvStart = nX+1+1;
        dxdvEnd   = nX+1+nX*nV;
        
        f       = m.f;
        dfdx    = m.dfdx;
        dfdv    = @dfdvSub;
        d2fdx2  = m.d2fdx2;
        d2fdvdx = @d2fdvdxSub;
        
        jDer = @der;
        jJac = @jac;
        if ~isempty(con.events)
            events  = con.events.events;
            jEvents = @eventFun;
        else
            jEvents = [];
        end
        
        % Derivative of [x; G; dxdv; dGdv] with respect to time
        function val = der(t, joint, u)
            u = u(t);            
            x = joint(1:nX);
            dxdv = reshape(joint(dxdvStart:dxdvEnd), nX,nV); % xv_ --> x_v
            
            % Compute derivative of dxdp
            dxdpdot = dfdx(t, x, u) * dxdv + dfdv(t, x, u); % f_x * x_v + f_v --> f_v
            
            % Sum continuous objective functions
            g = 0;
            dgdv = zeros(nV,1);
            for i = 1:nObj
                g = g + opts.ObjWeights(i) * obj(i).g(t, x, u);
                temp = vec(obj(i).dgdp(t, x, u)); % p_
                temp = [temp(:,opts.UseParams); sparse(nVX,1)]; % p_ --> v_   % remove rows for inactive parameters
                dgdv = dgdv + opts.ObjWeights(i) * (vec(obj(i).dgdx(t, x, u) * dxdv) + temp); % v_ + (_x * x_v --> _v --> v_) + v_ --> v_
            end
                        
            val = [f(t, x, u); g; vec(dxdpdot); dgdv];
        end
        
        % Jacobian of [x; G; dxdv; dGdv] derivative
        function val = jac(t, joint, u)
            u = u(t);            
            x = joint(1:nX); % x_
            dxdv = reshape(joint(dxdvStart:dxdvEnd), nX,nV); % x_v
            
            % Sum continuous objective gradients
            dgdx = zeros(1,nX); % _x
            d2gdxdv = sparse(nV,nX); % v_x
            for i = 1:nObj
                dgdx = dgdx + opts.ObjWeights(i) * obj(i).dgdx(t, x, u); % _x
                
                % Compute d/dx(dgdv)
                temp = obj(i).d2gdxdp(t,x,u); % Partial d2gdpdx p_x
                temp = [temp(opts.UseParams, :); sparse(nVX,nX)]; % p_x --> v_x
                d2gdxdv = d2gdxdv + opts.ObjWeights(i) * (dxdv.' * obj(i).d2gdx2(t,x,u) + temp); % v_x + v_x * x_x + v_x --> v_x
            end
            
            % Compute d/dx(dfdv)
            d2xdxdv = d2fdx2(t,x,u) * dxdv + d2fdvdx(t,x,u); % fx_v
            d2xdxdv = full(d2xdxdv); % fx_v
            d2xdxdv = reshape(d2xdxdv, nX,nX,nV); % f_x_v
            d2xdxdv = permute(d2xdxdv, [1,3,2]); % f_v_x
            d2xdxdv = reshape(d2xdxdv, nX*nV,nX); % fv_x
            
            % Combine
            val = [dfdx(t, x, u), sparse(nX,1+nX*nV+nV);
                            dgdx, sparse(1,1+nX*nV+nV);
                         d2xdxdv, sparse(nX*nV,1), kron(Ip, dfdx(t,x,u)), sparse(nX*nV,nV);
                         d2gdxdv, sparse(nV,1), kron(Ip, dgdx), sparse(nV,nV)];
        end
        
        % Modifies dfdp to relate only to the parameters of interest
        function val = dfdvSub(t, x, u)
            val = m.dfdp(t, x, u);
            val = [val(:, opts.UseParams) zeros(nX, nVX)];
        end
        
        % Modifies d2fdpdx to relate only to the parameters of interest
        function val = d2fdvdxSub(t, x, u)
            val = m.d2fdpdx(t,x,u);
            val = [val(:, opts.UseParams) zeros(nX*nX, nVX)];
        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