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: 
Here's a way to develop a series solution for a quasilinear ODE with analytic coefficients, $y^{(n)}=f(x,y,y',\dots,y^{(n-1)})$. The main idea to fold `Solve` over the series expansion of the ODE solve for the derivatives of $y$, something like this:

```
Fold[
  Flatten[{#1, Solve[#2 == 0 /. #1]}] &,
  ics,                                               (* initial conditions *)
  Series[ode /. Equal -> Subtract, {x, x0, n}][[3]]  (* series coefficients *)
  ]

```

Getting this idea to work for others without their having tweaking it takes some work. I've handled most of the things I've run into in the code below, but there are probably some things I have not tested. A malformed ODE may fail, probably ungracefully. Suggestions and comments are welcome if you run into other problems.

```
Clear[seriesDSolve];
seriesDSolve[ode_, y_, iter : {x_, x0_, n_}, ics_: {}] := 
 Module[{dorder, coeff},
  dorder = Max[0, Cases[ode, Derivative[m_][y][x] :> m, Infinity]];
  (coeff = Fold[
      Flatten[{#1, 
         Solve[#2 == 0 /. #1, 
          Union@Cases[#2 /. #1, Derivative[m_][y][x0] /; m >= dorder, 
            Infinity]]}] &,
      ics,
      Table[
       SeriesCoefficient[ode /. Equal -> Subtract, {x, x0, k}],
       {k, 0, Max[n - dorder, 0]}]
      ];
    Series[y[x], iter] /. coeff) /; dorder > 0
  ]

```

Example (constructed from a comment by the OP):

```
seriesDSolve[
 y''[x] + y'[x] == Sin[2 y[x]],
 y, {x, 0, 3},
 {y[0] -> a, y'[0] -> b}         (* symbolic initial values *)
 ]

```

![Mathematica graphics](https://i.stack.imgur.com/xuunF.png)

(Looking at it, I suppose it would be nice if the initial condition could be specified just like in `DSolve`. That would take some parsing and this approach evolved according to in-house needs -- that is to say, mine.)
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 *)

```