Question: ODE numeric integration

I need to integrate an ODE : I use dsolve.

As there is not an exact solution, I use numeric integration :

F := y''(x) + a(x)*y'(x) + b(x)*y(x) ;

EDF := {F = 0, y(t0) = y0, y'(t0) = dy0};

Snum := dsolve(EDF, y(x), type = numeric);

S  := x --> rhs(Snum(x)[2]) ;

I can now plot the graph of the solution S.

But suppose I need also the graph of the derivative S'.
As this function is solution of Y' + a*Y + b*S = 0,
I code

F1:=Y'(x) + a(x)*Y(x) + b(x)*S(x) ;

EDF1 := {F1 = 0, Y(t0) = dy0};

S1num := dsolve(EDF1, Y(x), type = numeric);

S1 := x --> rhs(S1num(x)[2]) ;

But Maple does not accept S(x) inside the définition of EDF1.
Thus, I must use a new function, say SS(x), which is an approximation of S(x) by interpolation.

Any suggestion ?

Thank you.