function sol = integrateSens(m, con, opts)

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

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

%TODO: allow for starting from steady-state

% Initial conditions [x0; vec(dxdv0)]
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
ic = [x0; vec([dxdp0, dxdx0])];

% 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, [], ic, con, opts.AbsTol(1:nX+nX*nV), opts.RelTol, stopTimes, 1, 1:nX);

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

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

        dxdvStart = nX+1;
        dxdvEnd   = nX+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; dxdv] 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
            
            % Combine
            val = [f(t, x, u); vec(dxdpdot)];
        end
        
        % Jacobian of [x; dxdv] derivative
        function val = jac(t, joint, u)
            u = u(t);            
            x = joint(1:nX); % x_
            dxdv = reshape(joint(dxdvStart:dxdvEnd), nX,nV); % x_v
            
            % 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,nX*nV);
                         d2xdxdv, kron(Ip, dfdx(t,x,u))];
        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