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: 
This method of developing a truncated solution can be done as below. I illustrate with an example that DSolve does not seem much to like.

```
ode = x''[t] - t*x'[t] + Sin[t] == 0;
initconds = {x'[0] == 1, x[0] == 0};

```

We create a differentail operator to create this ode.

```
odeOperator = D[#, {t, 2}] - t*D[#, t] + Sin[t] &;

```

Now set up our Taylor series as a symbolic expansion using derivatives of `x` evaluated at the origin. I use an order of 15 but that is something one would probably make as an argument to a function, if automating all this.

```
xx = Series[x[t], {t, 0, 15}];

```

Next apply the differential operator and add the initial conditions. Then find a solution that makes all powers of `t vanish.

```
soln = SolveAlways[Join[{odeOperator[xx] == 0}, initconds], t];

```

Let's look at the Taylor polynomial.

```
truncatedSol = Normal[xx /. soln[[1]]]

(* Out[500]= t + t^5/120 + t^7/1260 + (29 t^9)/362880 + (
 13 t^11)/1995840 + (2861 t^13)/6227020800 + (4649 t^15)/163459296000 *)

```

To assess how accurate it might be I will compare with `NDSolve`.

```
approxSol = NDSolve[Join[{ode}, initconds], x[t], {t, 0, 4}];

Plot[{truncatedSol, x[t] /. approxSol[[1]]}, {t, 0, 4}]

```

![enter image description here](https://i.stack.imgur.com/UsSAm.png)
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.)