Derivative of solution of differential equation with respect to parameter Find derivative with respect to $A$ of solution of differential equation
$$
\ddot x = {\dot x}^2 + x^3\tag1\label1
$$
with initial conditions $\{x(0)=0,\,\dot x(0)=A \}$ at $A=0$.
My attempt
We can find direct solution of $\eqref{1}$, but only implicitly (using substitution $p=\dot x^2$). Anyway, I prefer series:
$$
x=x_0 + Ax_1 + \ldots,\tag2\label2
$$
where $x_0, x_1, \ldots$ are functions of $t$. Now the question is posed as to find $x_1$. From initial conditions, $x_0(0)=0, x_1(0)=0, \dot x_0(0)=0, \dot x_1(0)=1$. Substituting $\eqref2$ into $\eqref1$, we have
$$
\ddot x_0 = \dot x_0^2 + x_0^3,\\
\ddot x_1 = 2\dot x_0\dot x_1 + 3x_0^2 x_1 \tag3
$$
But now I stuck. I can find $x_0$ only implicitly, so I cannot solve second equation for $x_1$. How to proceed from this ot how to use another methods?
 A: *

*OP is considering the 2nd order ODE 
$$\ddot{x} ~=~ \dot{x}^2 + x^3, \qquad x(0)~=~0,\qquad \dot{x}(0)~=~A.\tag{A}$$

*Notice that when $A=0$, we get the constant solution $$x_0~\equiv~ 0.\tag{B}$$

*OP's method. Expand the ODE (A) in a power series $$x~=~x_0 +Ax_1 + \frac{A^2}{2}x_2+ \ldots , \qquad x_n(0)~=~0,\qquad \dot{x}_n(0)~=~\delta_n^1,\tag{C}$$
in $A$. To first order in $A$, we get
$$\ddot{x}_1~=~0\qquad x_1(0)~=~0,\qquad \dot{x}_1(0)~=~1,\tag{D} $$
cf. OP's eq. (3b). The solution is 
$$ x_1(t)~=~t.\tag{E} $$

*Alternative method. A first integral to ODE (A) is$^1$
$$ \dot{x}^2~=~\left(A^2+\frac{3}{4}\right)e^{2x}-2V(x), \qquad x(0)~=~0,\tag{F}$$
where
$$ 2V(x)~:=~\frac{4x^3+6x^2+6x+3}{4}.\tag{G} $$ 

*Next expand the first integral (F) in a power series (C) To second order in $A$, we get 
$$ \dot{x}_1^2~=~1, \qquad  x_1(0)~=~0, \tag{H}$$
with solution
$$ x_1(t)~=~\pm t.\tag{I} $$
The minus sign branch in eq. (I) is discarded by comparing with initial conditions.
--
$^1$ In fact, eq. (A) is the Euler-Lagrange eq. for the Lagrangian
$$ L(x,\dot{x})~=~e^{-2x}\left(\frac{1}{2}\dot{x}^2-V(x)\right). \tag{J}$$
Since there is no explicit time dependence, the corresponding energy function
$$h(x,\dot{x})~:=~\dot{x}\frac{\partial L}{\partial \dot{x}}-L
~=~e^{-2x}\left(\frac{1}{2}\dot{x}^2+V(x)\right)\tag{K}$$
is conserved, cf. Noether's theorem.
