How to estimate settling time of an overdamped system?

Well, let's solve this in a more general case. We have the following transfer function (assuming real positive value for $epsilon$):

$$mathcal{H}left(text{s}right):=frac{text{Y}left(text{s}right)}{text{X}left(text{s}right)}=frac{1}{left(text{s}+epsilonright)left(text{s}+2epsilonright)}tag1$$

When we look at the step-response we are using $text{X}left(text{s}right)=mathcal{L}_tleft[thetaleft(tright)right]_{left(text{s}right)}=frac{1}{text{s}}$, so the output is given by:

$$text{Y}left(text{s}right)=frac{1}{text{s}}cdotfrac{1}{left(text{s}+epsilonright)left(text{s}+2epsilonright)}tag2$$

Using inverse Laplace transform, we find:

$$text{y}left(tright)=mathcal{L}_text{s}^{-1}left[frac{1}{text{s}}cdotfrac{1}{left(text{s}+epsilonright)left(text{s}+2epsilonright)}right]_{left(tright)}=frac{expleft(-2epsilon tright)left(expleft(epsilon tright)-1right)^2}{2epsilon^2}tag3$$

It is not hard to show that when $ttoinfty$ (assuming real positive value for $epsilon$), we get:

$$lim_{ttoinfty}text{y}left(tright)=frac{1}{2epsilon^2}tag4$$

Now, for the settling time, we want to find the time $t$ when $text{n}text{%}$ of the final value is reached:

$$text{y}left(t_text{n}right)=frac{text{n}text{%}}{100}cdotfrac{1}{2epsilon^2}Longleftrightarrowspace t_text{n}=dotstag5$$

Solving that gives:

$$t_text{n}=frac{1}{epsilon}cdotlnleft(frac{10}{10-sqrt{text{n}text{%}}}right)tag6$$

For systems with real left-half plane poles, you can usually estimate it by only considering the dominant pole (the pole with the lowest frequency). In your case this would be $p_d=-2$. The result gets more accurate as the non-dominant pole ($p_{nd}$) moves further away from the dominant pole.

By only considering the dominant pole, you get a rather simple equation:

$$begin{align} frac{5}{4}cdot e^{-2t}&=0.02cdot frac{5}{8} \ t &= -frac{1}{2}cdot lnleft( 0.02cdot frac{5}{8}cdot frac{4}{5} right) approx 2.30258509299s\ end{align}$$

The idea is that the non-dominant pole at $p_{nd}=-4$ leads to a term $e^{-4t}$ which will be damping so quickly that it doesn't affect the overall settling time. The advantage is the simplicity of the equation, and the fact that is actually a pretty common occurrence to have a very dominant pole and far-away non-dominant poles in electronic circuits.

In your specific case, it is possible to analytically calculate the settling time. The time it takes for the time-dependent terms to damp to 2% of the final value can be calculated using (similar to Andy's answer, but using the absolute value):

$$begin{align} left| e^{-4t}-2cdot e^{-2t} right| &=0.02 \ &Updownarrow (y=e^{-2t})\ y^2-2cdot y &= pm 0.02 \ &Updownarrow (text{There are 4 distinct solutions, but I only take the relevant one}) \ y = e^{-2t} &= 1 - frac{7}{5sqrt{2}} \ &Updownarrow \ t &= -frac{1}{2}cdot lnleft(1 - frac{7}{5sqrt{2}}right) approx 2.30006613189s end{align}$$

So a factor of 2 for $p_{nd}/p_d$ leads to about an error of 0.1% on the calculated settling time when using the dominant-pole approximation. Whether or not this is sufficient I leave to you.

Yes, your inverse Laplace calculation is correct.

The final steady state value will be 5/8 - this is the DC value after a long length of time. So, you are really looking for the rest of the equation to fall in magnitude to 2% of 5/8: -

$$dfrac{5}{8}e^{-4t} - dfrac{5}{4}e^{-2t} = dfrac{5}{8}cdot text{0.02}$$

$$=dfrac{8}{8}e^{-4t} - dfrac{8}{4}e^{-2t} = dfrac{8}{8}cdot text{0.02}$$

$$= e^{-4t} - 2e^{-2t} = 0.02$$

Does that help?