Question: How do I find a series solution to an ODE? 
I do not mean taking the Taylor series of an exact solution; I want to solve nasty nonlinear differential equations locally via plug and chug. Surely, that is a built-in option in *Mathematica*.

Answer: 
[J.M.](https://mathematica.stackexchange.com/users/50) in a comment above alluded to this possibility. Using [`Series`](http://reference.wolfram.com/mathematica/ref/Series.html) and [`DifferentialRoot`](http://reference.wolfram.com/mathematica/ref/DifferentialRoot.html) seems like a great way to expand linear differential equations in series. It seems worth putting down in an answer due to its simplicity:

```
lde = {y''[x] - x*y'[x] + Sin[x] == 0, y'[0] == 1, y[0] == 0};
Series[DifferentialRoot[Function @@ {{y, x}, lde}][x], {x, 0, 15}]

(* x + x^5/120 + x^7/1260 + (29 x^9)/362880 + (13 x^11)/1995840 +
    (2861 x^13)/6227020800 + (4649 x^15)/163459296000 + O[x]^16 *)

```
On Mathematica `11.3` using a new function `AsymptoticDSolveValue`

Linear ODE:

```
AsymptoticDSolveValue[{y''[x] - x*y'[x] + Sin[x] == 0, y'[0] == 1, 
y[0] == 0}, y[x], {x, 0, 10}]

(* x + x^5/120 + x^7/1260 + (29 x^9)/362880  *)

```

Quasilinear ODE:

```
AsymptoticDSolveValue[{y''[x] + y'[x] == Sin[2*y[x]], y[0] == a, 
y'[0] == b}, y[x], {x, 0, 5}]

(* a + b x + 1/6 x^3 (b + 2 b Cos[2 a] - Sin[2 a]) + 
1/2 x^2 (-b + Sin[2 a]) + 
1/24 x^4 (-b - 4 b Cos[2 a] + Sin[2 a] - 4 b^2 Sin[2 a] + 
2 Cos[2 a] Sin[2 a]) + 
1/120 x^5 (b + 6 b Cos[2 a] - 8 b^3 Cos[2 a] + 4 b Cos[2 a]^2 - 
Sin[2 a] + 16 b^2 Sin[2 a] - 4 Cos[2 a] Sin[2 a] - 
12 b Sin[2 a]^2)*)

```

Nonlinear ODE:

```
AsymptoticDSolveValue[{y'[x] == 1 + y[x]^2, y[0] == 0}, y[x], {x, 0, 10}]

(* x + x^3/3 + (2 x^5)/15 + (17 x^7)/315 + (62 x^9)/2835 *)

```