A natural approach is to start from the special case where $A(t)=0$ for every $t\geqslant0$. Then the ODE reads y'(t)=B(t) hence y'(t)\geqslant0 for every $t\geqslant0$ hence $y$ is nondecreasing. Since $y(0)\geqslant0$, this proves that $y(t)\geqslant0$ for every $t\geqslant0$.
One can deduce the general case from the special one. This leads to consider $z(t)=C(t)y(t)$ and to hope that z'(t) is a multiple of the LHS of the ODE. As you know, defining $C(t)=\exp\left(\int\limits_0^tA(s)\mathrm ds\right)$ fits the bill since (C\cdot y)'=C\cdot (y'+A\cdot y) hence one is left with the ODE z'=C\cdot B.
Since $C(t)>0$, $C(t)\cdot B(t)\geqslant0$ hence the special case for the RHS $C\cdot B$ yields the inequality $z(t)\geqslant z(0)$ and since $z(0)=y(0)\geqslant0$, one gets $y(t)\geqslant C(t)^{-1}y(0)\geqslant0$ and the proof is finished.
Edit The same trick applies to more general functions $A$ and $B$. For example, replacing $A(t)$ by $A(t,y(t))$ and $B(t)$ by $B(t,y(t))$ and assuming that $B(s,x)\geqslant0$ for every $s\geqslant0$ and every $x\geqslant0$, the same reasoning yields $y(t)\geqslant C(t)^{-1}y(0)\geqslant0$ with $C(t)=\exp\left(\int\limits_0^tA(s,y(s))\mathrm ds\right)$.
The only modification is that one must now assume that $t\geqslant0$ belongs to the maximal interval $[0,t^*)$ where $y$ may be defined. Solutions of such differential equations when $A$ does depend on $y(t)$ may explode in finite time hence $t^*$ may be finite but $y(t)\geqslant0$ on $0\leqslant t.