Overshoot computations not matching
Computing the partial fraction expansion,

$$H (s) = \frac{2.302 s + 0.3548}{s^{3} + 0.739 s^{2} + 3.223 s + 0.3548} \approx \frac{2.2861 - 0.0309306 s}{s^2 + 0.626454 s + 3.1525} + \frac{0.0309306}{s + 0.112546}$$

For the time being, let us neglect the step response of the 1st order subsystem. Note that the 2nd order subsystem has the following form

$$\pm \gamma (s - z) \left(\frac{\omega_n^2}{s^2 + 2 \zeta \omega_n s + \omega_n^2}\right)$$

where $\gamma \in \mathbb R$ is the gain (or attenuation) and $z \in \mathbb R$ is a (finite) zero. However, the formula you used to calculate the step response overshoot is only applicable for 2nd order systems of the form

$$G (s) := \frac{\omega_n^2}{s^2 + 2 \zeta \omega_n s + \omega_n^2}$$

whose zeros are at infinity. To summarize, the overshoot is smaller than you expected because:


you neglected the 1st order subsystem.
you neglected the attenuation $\gamma \neq 1$.
you did not subtract $\gamma$ times the derivative of the step response of $G (s)$ at the peak time.


In any case, not neglecting the above would only provide one estimate of the overshoot. To compute a satisfactory approximation of its exact value, you can take the inverse Laplace transform, differentiate and find where the derivative vanishes.