I have discussed in another post how to work with state-space models to simulate the response of linear systems using Matlab’s lsim. It is however also possible to approach this problem from a transfer-function point of view without having to mess with state-space at all.
In order to appreciate what is going on one needs to look at systems of order higher than one. How about the generic second-order system modeled by the ODE
$$\ddot{y}(t) + a_1 \dot{y}(t) + a_2 y(t) = b_2 u(t)$$
The transfer-function $G$ in this case is
$$G(s) = \frac{b_2}{s^2 + a_1 s + a_2}$$
which is easily obtained with the help of the Laplace transform (See Chapter 3). More complex relationships involving derivatives of $u$ can be handled similarly as long as the resulting transfer-function remains proper (See Chapter 5).
The Laplace transform also provides a way to formally incorporate initial conditions. Using the differentiation in time property from Table 3.2 one can write
$$\mathcal{L} \{ y^{(n)}(t) \} = s^n F(s) – s^{n-1}y(0) – \cdots – y^{(n-1)}(0)$$
which applied to the above transfer-function provides
$$s^2 Y(s) – s y(0) – \dot{y}(0) + a_1 (s Y(s) – y(0)) + a_2 Y(s) = b_2 U(s)$$
Rearrange
$$(s^2 + a_1 s + a_2) Y(s) = b_2 U(s) + (s + a_1) y(0) + \dot{y}(0)$$
and solve for $Y(s)$ to obtain
$$Y(s) = G(s) U(s) + G_0(s)$$
where $G$ is the transfer-function you are familiar with and
$$G_0(s) = \frac{(s + a_1) y(0) + \dot{y}(0)}{s^2 + a_1 s + a_2}$$
is a rational function related to the initial conditions, in this case $y(0)$ and $\dot{y}(0)$. At this point one could simply invert $G_0$ and use linearity to add to the response. This can be done symbolically or by realizing that
$$\mathcal{L}^{-1} \{ G_0(s) \} = \mathcal{L}^{-1} \{ G_0(s) \times 1 \} = G_0 \delta(t)$$
so that this component can be calculated as an impulse response, that is
$$y(t) = G u(t) + G_0 \delta(t)$$
Because $\delta(t)$ is not a regular function you cannot use lsim to calculate that component of the response. Instead Matlab has the specialized function impulse that you can use. For example
y = lsim(G, u, t) + impulse(G0, t)
will do the job. Note that it is important to provide the vector t to the function impulse otherwise it will not return a vector compatible with the first vector calculated by lsim.
These ideas can be extended to higher-order systems and a formula for the numerator of $G_0(s)$ can be constructed if wanted. Such formula is implicitly related to one state-space realization and the procedure outlined here.