Here is the damped oscillator and its generic solution:

de1 = x''[t]+ 2 γ x'[t] + ω0^2 x[t] == 0
s1=DSolve[de1,x[t],t]

A typical problem might involve determining C[1] and C[2] from some initial value data, in which case we are interested in the position and velocity at time t=0:

x[t]/. s1[[1]]  /. t->0
D[x[t]/. s1[[1]],t]/.t->0

DSolve does the Solve[]ing for us:

s1=(DSolve[{de1,x[0]==x0,x'[0]==v0},x[t],t]//FullSimplify)[[1]]

with some choice of numerical parameters:

p1= {ω0->1,γ->1/10,x0->1,v0->0}

x1[t_] = x[t]/.s1/.p1

ComplexExpand[x1[t]]// Simplify[#,Assumptions->{t>0}]&

x1[t_] = %

Express the period in terms of givens, and plot one cycle:

T = 2 π / Sqrt[ω0^2-γ^2]/.p1
Plot[x1[t],{t,0,T}]

Note the zero comes “late” (i.e. more than π/2 in the cycle):

FindRoot[x1[t],{t,1}]
(t/.%)/(π/2)

Here is the energy:

e[t_] = (1/2) (x1'[t])^2 + (1/2)(ω0/.p1)^2 x1[t]^2//Simplify

LogPlot[e[t],{t,0,50}]

Now add a driving force

de2 = x''[t]+ 2 γ x'[t] + ω0^2 x[t] == Cos[ω t]

The generic solution looks like:

s=DSolve[de2,x[t],t]//Simplify

Let’s pick out that long time non-decaying piece:

x2[t_] = x[t]/. s[[1]] /. {C[1]->0,C[2]->0}

Here is the straightforward but painful way of finding the amplitude of oscillations:

peaktime=Solve[x2'[t]==0,t];
amp=x2[t]/.peaktime[[1]]//FullSimplify[#,Assumptions->{γ>0,ω>0,ω0>0}]&

Ouch!  A Mathematica bug?!  The correct answer is:

amp= ((+4 (γ^2) (ω^2)+((((ω^2)-((ω0)^2)))^2))/(((4 (γ^2) (ω^2)+((((ω^2)-((ω0)^2)))^2)))^(3/2)))

Plot[amp/.{ω0->1,γ->1/10},{ω,0,3}]

Pick some driving frequency and solve:

p2= {ω0->1,γ->1/10,ω->3,x0->1,v0->1};
s2=DSolve[{de2 , x[0]==x0,x'[0]==v0}/.p2,x[t],t][[1]]

x2[t_] = Simplify[x[t] /. s2 ]

Plot[x2[t],{t,0,70}]

Plot[x2[t],{t,30,100}]

Here we include the driving force in the plot, scaled to fit:

Plot[{x2[t],.3*Cos[ω t]/.p2},{t,0,20}]

Plot[{x2[t],.2*Cos[3 t]},{t,90,94}]

Find the phase shift:

Manipulate[Plot[{x2[t],.2*Cos[(ω/.p2)t+ϕ]},{t,90,94}],{ϕ,0,2π}]