Perturbation series for projectile problem

\(TraditionalForm
    \`Generation\ of\ a\ regular\ perturbation\ expansion\ by\ series\ 
      expansion\ for\ the\ problem\)\(TraditionalForm
    \`\(\(d\^2\) y\)\/\(d\ t\^2\) = 
      \(-\ \(1\/\((1 + ϵ\ y)\)\^2\)\)\ \ , \ \ y[0] = 0, \ \ 
    \(y'\)\ [0] = 1.\)

Determining the first three terms of the perturbation series expansion

We can start by assuming an expansion of the following form (here first three terms only - we will extended this below to any order)

y[t_] := (Subscript[y, 0])[t] + ϵ (Subscript[y, 1])[t] + ϵ^2 (Subscript[y, 2])[t]

After shifting all terms to one side, the differential equation becomes

de  = y''[t] + 1/(1 + ϵ y[t])^2

We Taylor expand this in ϵ

ser = Series[de, {ϵ, 0, 2}]

The system of ODEs to solve becomes

odes = LogicalExpand[ser == 0]

It is safest to store the ODEs in a list rather than as a logical statement when we need to refer to them in turn

ode = Table[odes[[i]], {i, 1, 3}]

Now solve in turn for (Subscript[y, 0])[t], (Subscript[y, 1])[t], and (Subscript[y, 2])[t]. Note the boundary conditions differ between the first- and subsequent cases. To avoid infinite recursions, we need to use an intermediate variable  temp  when defining the functions (Subscript[y, k])[t_] from the solutions of the ODEs.

sol = 
      DSolve[{ode[[1]], (Subscript[y, 0])[0] == 0, (Subscript[y, 0])'[0] == 1}, (Subscript[y, 0])[t], t]; 
    temp = (Subscript[y, 0])[t] /. sol[[1]];  (Subscript[y, 0])[t_] = temp

sol = 
      DSolve[{ode[[2]], (Subscript[y, 1])[0] == 0, (Subscript[y, 1])'[0] == 0}, (Subscript[y, 1])[t], t]; 
    temp = (Subscript[y, 1])[t] /. sol[[1]];  (Subscript[y, 1])[t_] = temp

sol = 
      DSolve[{ode[[3]], (Subscript[y, 2])[0] == 0, (Subscript[y, 2])'[0] == 0}, (Subscript[y, 2])[t], t]; 
    temp = (Subscript[y, 2])[t] /. sol[[1]];  (Subscript[y, 2])[t_] = temp

\(We\ have\ now\ arrived\ at\ \ y[t, ϵ] = 
            t - t\^2\/2 + \((t\^3\/3 - t\^4\/12)\)\ ϵ + 
              \((\(-\(t\^4\/4\)\) + \(11\ t\^5\)\/60 - \(11\ t\^6\)\/360)\)\ 
                ϵ\^2 + O[ϵ]^3.

          
          This\ particular\ code\ becomes\ somewhat\ tedious\ when\ we\ want\ 
            many\ terms.\ Below\ is\ a\ version\ that\ \)automatically loops to give any number of terms.

Code for unlimited number of terms

Remove["Global`*"]

The code below is essentially identical to the one above - instead of solving for the first three expansion functions one-by-one, we loop over an arbitrary number of terms. Suppose we want to get all terms up to (and including)  ϵ^n, in case of n = 6:

n = 6

\(TraditionalForm
    \`This\ means\ determining\ all\ the\ functions\ y\_k[t]\ , \ k\  = \ 0, 
    1, ..., n\ defined\ through\)

y[t_] := Sum[(Subscript[y, k])[t] ϵ^k , {k, 0, n}]

The differential equation is again  de = 0  where

de  = y''[t] + 1/(1 + ϵ y[t])^2

We Taylor expand this in ϵ

ser = Series[de, {ϵ, 0, n}]

and get the system of ODEs through

odes = LogicalExpand[ser == 0]; 
    ode = Table[odes[[i + 1]], {i, 0, n}]

Now, we loop over the successive ODEs.

\(Do\ [
      sol = \ DSolve[{ode\ [\([k + 1]\)]	, y\_k[0] == 0, 
            \(y\_k'\)[0] == If[k == 0, 1, 0]}, y\_k[t], t]\ ; 
	
      temp = y\_k[t] /. sol[\([1]\)]; \ y\_k[t_] = temp, 

      	\ \ \ \ {k, 0, n}]\)

The expansion up to and including  (ϵ^n)  becomes

y[t]