function sol = integrateSys(m, con, opts)
%% Run integration for each experiment
nX = m.nX;

% Construct system
[xDer, xJac, xEvents] = constructSystem();

%TODO: allow for starting from steady-state

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

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

% Integrate x over time
sol = accumulateSol(xDer, xJac, xEvents, [], ic, con, opts.AbsTol(1:nX), opts.RelTol, stopTimes, 1, 1:nX);

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% The system for integrating x %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [xDer, xJac, xEvents] = constructSystem()
        f    = m.f;
        dfdx = m.dfdx;
        
        xDer = @der;
        xJac = @jac;
        
        if ~isempty(con.events)
            events  = con.events.events;
            xEvents = @eventFun;
        else
            xEvents = [];
        end
        
        % Derivative of x with respect to time
        function val = der(t, x, u)
            u   = u(t);
            val = f(t, x, u);
        end
        
        % Jacobian of x derivative
        function val = jac(t, x, u)
            u   = u(t);
            val = dfdx(t, x, u);
        end
        
        % Event function for system
        function [val, dir, term] = eventFun(t, x, u)
            u = u(t);
            [val, dir, term] = events(t, x, u);
        end
    end
end