t=0:0.01:150*pi;
omp=0.9;  % frequency of damped oscillator
omf=0.9;  % frequency of external forcing
Q=20;     % Quality factor
gam=omp/Q; % dissipation factor
A=1.5;X=1; % steady state and transient amplitudes 
del=atan(omf*gam/(omp^2-omf^2)); 
alpha=1*pi/2;  % initial phase of free oscillation

x=A.*cos(omf*t-del)+X.*exp(-gam/2*t).*cos(omp*t+alpha);
plot(t,x)
axis([0 max(t) -2 2])