We can find many step-response-based identification methods in the literature for a non-oscillatory system (a system with solid damping or an overdamped system). One of them is Harriott's method, which can be considered a graphical method (some call it a semi-graphical method). The following lecture note is an excellent reference to Harriott's method.
The following lecture note is an excellent reference to Harriott’s method.
The concept is quite difficult to understand. Mainly, Harriott introduced a model with two time constant $\tau_1$ and $\tau_2$, where $\tau_1>\tau_2$. The model is as follows:
$$ G(s)=\frac{K_{ss}}{(\tau_1 s+1)(\tau_2 s + 1)}e^{-\theta s} $$
After some careful investigations, Harriott found that $t_{72}=\theta+1.25(\tau_1+\tau_2)$. Next, based on an empirical plot generated from the step responses of $t=0.4t_{72}+0.6\theta$, Harriott sould find the relationship between $\tau_1$ and $\tau_2$.
Let us for now ignore the details on the Harriott’s method and directly translate the lecture note above into MATLAB codes.

There are at least six other methods if we refer to this paper. This one is probably the most recent. They are all interesting from an industrial perspective since the implementations are simple and quick even though the results are most likely to be less accurate than the optimization-technique-based results.