How to linearize an expression automatically?

I'd use Series :

f[a1_, a2_] = a1^2 a2;

(Series[f[q1 + \[Epsilon] dq1, q2 + \[Epsilon] dq2], {\[Epsilon], 0, 1}] // Normal) 
  /. \[Epsilon] -> 1
(* dq2 q1^2 + 2 dq1 q1 q2 + q1^2 q2 *)