There is no builtin function I could find in Mathematica to do this. But one way is the following. We will have a function that returns the coefficients of the ode and the RHS. i.e. given an ode \(A y''+ By' + C y=f(x)\), it will return a list \(\{A,B,C,f\}\). Next to check if the ode is linear or not, all what we have to do is check that \(A,B,C\) do not have \(y\) in them. This include \(y,y',y''\) and so on.
So for a first order ode, the result of calling this function, which we will get getODEcoeff is always a list. The Length of the list is 2 more than the order of the ode. If a coefficient is missing, zero is placed in that slot. For example, for 3 y''[x] + x y[x]==Sin[x] then getODEcoeff[ode,y[x],x] will return {3,0,x,Sin[x] and for x y''''[x] + x y[x]==Sin[x] it will return {x,0,0,0,x,Sin[x]} and for first order
ode y'[x]+f[x]*y[x]==9 it will return {1,f[x],9}.
Now we simply check that all element of this list have no \(y\) in them. Here is the function
(to do)