dsolve(diff(diff(y(x),x),x) = -(a*p*x^b+q)/x/(a*x^b-1)*diff(y(x),x)-(a*r*x^b+s)/x^2/(a*x^b-1)*y(x),y(x))