I suppose that ParametricNDSolve could present another approach.
func = ParametricNDSolveValue[{a'[t] == w[t], w'[t] == -k w[t] - Sin[a[t]] + l Cos[q[t]], q'[t] == wd,
a[0] == a0, w[0] == w0, q[0] == 0}, {a, w, q}, {t, 0, 100}, {a0, w0, k, l, wd}]
pend[a0_, w0_, k_, l_, wd_] := func[a0, w0, l, k, wd]
Plot[#[t] & /@ pend[1., 2., 3., 4., 5.][[1 ;; 2]], {t, 0, 100}]
Cheers,
Marco