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)
[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 *)

```