# Plotting a complete response to an oscillation; already have particular solution, not sure how to solve for constants in homogeneous solution

I am writing a code that plots an oscillator's response to a periodic excitation. The response consists of the following:

x(t) = xh + xp

where xh is the homogeneous solution and xp is the particular solution. I have already successfully figured out how to plot the particular solution, which is represented as a Fourier series with coefficients a0, an, and bn.

I know that the homogeneous solution is of the form:

but I am not sure how to solve for the constants. The initial conditions are that x[0]=0.2 and x'[0]=-0.05 (this is applied for x(t) = xh + xp). I only really know how to solve for these constants when given a differential equation since Mathematica has a function for that, but I am not sure how to ask it to solve for c1 and c2 with both the Fourier series added to the pasted solution above. Everything in the code up to where I define xh correctly works (I apologize for the greek characters in this code I couldn't figure out how to make them look right here, this is my first time posting):

ClearAll["Global*"]
m = 15;
k = 960;
c = 24;
ωn = Sqrt[k/m]
ζ = c/(2*m*ωn)
ωd = ωn*Sqrt[1 - ζ^2]
Te = 4;
F = 60;
F0 = ((F*Te)/4)*t;
F1 = 90 - 30*t;
F2 = 0;
ω = (2*Pi)/Te
r = ω/ωn
MF = 1/((1 - (n*r)^2)^2 + (2*ζ*n*r)^2)^(1/2)
θ = ArcTan[1 - (n*r)^2, 2*ζ*n*r];
a0 = (2/Te)*(Integrate[F0, {t, 0, Te/4}] + Integrate[F1, {t, Te/4, (3*Te)/4}] + Integrate[F2, {t, (3*Te)/4, Te}])
an = (2/Te)*(Integrate[F0*Cos[(2*n*Pi*t)/Te], {t, 0, Te/4}] + Integrate[F1*Cos[(2*n*Pi*t)/Te], {t, Te/4, (3*Te)/4}] + Integrate[F2*Cos[(2*n*Pi*t)/Te], {t, (3*Te)/4, Te}])
a[n_] := Simplify[an, Element[n, Integers]]
a[n]
bn = (2/Te)*(Integrate[F0*Sin[(2*n*Pi*t)/Te], {t, 0, Te/4}] + Integrate[F1*Sin[(2*n*Pi*t)/Te], {t, Te/4, (3*Te)/4}] + Integrate[F2*Sin[(2*n*Pi*t)/Te], {t, (3*Te)/4, Te}])
b[n_] := Simplify[bn, Element[n, Integers]]
b[n]
xss[t_] := a0/(2*k) + Sum[(a[n]/k)*MF*Cos[n*ω*t - θ] + (b[n]/k)*Sin[n*ω*t - θ], {n, 1, 100}]
xh[t_] = Exp[(-ζ)*ωn*t]*(c1*Cos[ωd*t] + c2*Sin[ωd*t])
x[t_] = xss[t_] + xh[t_]
Solve[x[0] == 0.2, Derivative[1][x][0] == -0.05, {c1, c2}]
Plot[x[t], {t, 0, 12}]


In general you would use Solve but here it is straightforward:

C1=0.2-xp[0]

and

I let you write the initial equation for C2.

ClearAll["Global*"]
m = 15;
k = 960;
c = 24;
ωn = Sqrt[k/m];
ζ = c/(2*m*ωn);
ωd = ωn*Sqrt[1 - ζ^2];
Te = 4;
F = 60;
F0 = ((F*Te)/4)*t;
F1 = 90 - 30*t;
F2 = 0;
ω = (2*Pi)/Te;
r = ω/ωn;
MF = 1/((1 - (n*r)^2)^2 + (2*ζ*n*r)^2)^(1/2);
θ = ArcTan[1 - (n*r)^2, 2*ζ*n*r];
a0 = (2/Te)*(Integrate[F0, {t, 0, Te/4}] +
Integrate[F1, {t, Te/4, (3*Te)/4}] +
Integrate[F2, {t, (3*Te)/4, Te}]);
an = (2/Te)*(Integrate[F0*Cos[(2*n*Pi*t)/Te], {t, 0, Te/4}] +
Integrate[F1*Cos[(2*n*Pi*t)/Te], {t, Te/4, (3*Te)/4}] +
Integrate[F2*Cos[(2*n*Pi*t)/Te], {t, (3*Te)/4, Te}]);
a[n_] := Simplify[an, Element[n, Integers]];
bn = (2/Te)*(Integrate[F0*Sin[(2*n*Pi*t)/Te], {t, 0, Te/4}] +
Integrate[F1*Sin[(2*n*Pi*t)/Te], {t, Te/4, (3*Te)/4}] +
Integrate[F2*Sin[(2*n*Pi*t)/Te], {t, (3*Te)/4, Te}]);
b[n_] = Simplify[bn, Element[n, Integers]];
xss[t_] :=
a0/(2*k) +
Sum[(a[n]/k)*MF*Cos[n*ω*t - θ] + (b[n]/k)*
Sin[n*ω*t - θ], {n, 1, 5}];
xh[t_] :=
Exp[(-ζ)*ωn*t]*(c1*Cos[ωd*t] +
c2*Sin[ωd*t]);
x[t_] := xss[t] + xh[t];
sol = NSolve[{x[0] == 0.2, (D[x[t], t] /. t -> 0) == -0.05}, {c1, c2}];
Plot[x[t] /. sol, {t, 0, 12}]