(* Rieder & Busby, 1986, Exa. 5.3, p 109 --- 2011-Jun MCasquilho *)
a = 3; alpha = 1.2; b = 1; ti = 0; tf = 0.5; y0 = 1;
sol = NDSolve[{y'[t] == a y[t]^alpha + b t, y[0] == y0}, y, {t, ti, tf}];
Show[Plot[y[u] /. sol, {u, ti, tf}, PlotStyle -> Automatic],
     Plot[y'[u] /. sol, {u, ti, tf}, PlotStyle -> Automatic] ]
f = sol[[1,1,2]]; f[tf]

(* Plot[Evaluate[{y[x],y'[x]} /. sol], {x, ti, tf}, PlotStyle -> Automatic] *)