You need specific boundary conditions in order to plot, and specific w_n value.
syms x_d(t_d)
syms w_n
ode = diff(x_d,t_d,2)+diff(x_d,t_d)-x_d == cos(w_n*t_d);
x_dSol(t_d) = dsolve(ode)
x0 = randi([-9 9]), dx0 = randi([-9 9]), w_n0 = randi([1 300])
dx_d = diff(x_d,t_d);
ic = [x_d(0) == x0, dx_d(0) == dx0]
ode1 = subs(ode, w_n, w_n0)
%you SHOULD be able to solve with two initial conditions, but dsolve does not
%like it. So solve with x(0) boundary condition, then differentiate
%and solve for the constant
x_dSol2 = dsolve(ode1, ic(1))
const_var = setdiff(symvar(x_dSol2), t_d);
const_sol = solve(subs(diff(x_dSol2, t_d), t_d, 0), const_var)
x_dSol3 = simplify(subs(x_dSol2, const_var, const_sol))
fplot(x_dSol3, [0 10])
