from pylab import *
from numpy import *

m = 1.
k = 1.
gamma = 0.1
F = 1
Omega = 1.

x0 = 1.
v0 = 0.

beta = gamma/(2*m)
omega = sqrt(k/m - beta**2)

A=F/m/((k/m-Omega**2)**2 + (gamma*Omega/m)**2)
if(k/m-Omega**2 == 0.):
  Fi = pi/2
else:
  Fi=arctan(gamma*Omega/m/(k/m-Omega**2))

fi=arctan(beta/omega +(A*Omega*sin(Fi)+v0)/(A*cos(Fi)-x0)/omega)
a = (-A*cos(Fi)+x0)/cos(fi)

t = linspace(0,100,10000)
x = A*cos(Omega*t + Fi)+a*exp(-beta*t)*cos(omega*t + fi)
    
plot(t,x)