An interesting ODE I've been thinking about approaches to solve
$$\frac{\textrm{d}^{3}y}{\textrm{d}x^{3}}=\textrm{e}^{-y(x)}.$$
My initial thought was to set $y(x)=\ln(z(x))$ in order to obtain an ODE relating $z$ to $x$. However, this in fact makes things much worse (or at least much messier).
Any insights much appreciated.
Thanks
 A: Your equation is autonomous. The standard method is to let $u(y)=y'(x)$, then carefully writing out the equation we have
$$
u(u')^2+u^2 u''=e^{-y}
$$
Now $y$ is the independent variable. A particular solution is $u_p=\left( \frac{9}{2} \right)^{1/3}e^{-y/3}$, letting $u=u_p+v$ just produces another awful nonlinear equation for $v$. At least we have one solution to the original equation: $y_p(x)=3 \ln(x-c_1)+\ln(6)$. In the absence of some clever integrating factor trick, we can look at perturbative solutions
$$
y'''=\varepsilon e^{- y}
$$
When $\varepsilon=0$, we have $y_0(x)=Ax^2+Bx+C$. With the ansatz
$$
y(x)=\sum\limits_{n=0}^\infty y_n(x)\varepsilon^n
$$
Choosing $y(0)=1, \ y'(0)=-1, \ y''(0)=0$, we have $y_0=1-x$, and we can find the next term by matching powers of $\varepsilon^\dagger$
$$
y_1'''=e^{-y_0}=e^{-1+x}
$$
Which is easily integrated
$$
y_1(x)=-\frac{1}{2e}(2-2e^x +2x+x^2)
$$
Here is a plot of the numerical integration and the two term series $y_0+y_1$:

$\dagger$ To get the higher order terms, you'd want a systematic way of getting down the right coefficients of $\varepsilon$ from $\sum_n y_n \varepsilon^n$ that's in the exponential. For example, at the next order, when equating $\varepsilon^2$ terms, we have on the right
$$
\varepsilon e^{-y_0-\varepsilon y_1 - \dots}=\varepsilon e^{-y_0} e^{-\varepsilon y_1} \dots =\varepsilon e^{-y_0}(1-\varepsilon y_1) \dots
$$
So that the $\varepsilon^2$ coefficient on the right is $-y_1 e^{-y_0}$, and the next term in our perturbation series is given by solving
$$
y_2'''=-y_1 e^{-y_0}
$$
Update: while I do not have a closed form for the coefficients of $\varepsilon^n$, they can be found diagrammatically, as they are in correspondence with the integer partitions of $n$, or equivalently the Young diagrams of $n$ blocks. Consider
$$
e^{-\sum \varepsilon^n y_n}=e^{-y_0} e^{-\varepsilon y_1} e^{-\varepsilon^2 y_2} e^{-\varepsilon^3 y_3}  \dots \\
=e^{-y_0}\left(1-\varepsilon y_1 + \frac{\varepsilon^2 y_1^2}{2} \dots \right) \left(1-\varepsilon^2 y_2 + \frac{\varepsilon^4 y_2^2}{2} \dots \right)\left(1-\varepsilon^3 y_3 + \frac{\varepsilon^6 y_3^2}{2} \dots \right)\dots
$$
The coefficient of $\varepsilon^n$ (from the exponential, excluding the single factor outside the exponential) will be given by all the ways in which expanding the parentheses out can produce $\varepsilon^n$. For example, the coefficients of $\varepsilon^4$ would be represented by the diagrams of four blocks:

Where I'm letting a column of length $n$ = one factor of $y_n$. Thus the fifth order equation reads
$$
y_5'''=e^{-y_0} \left( \frac{(-y_1)^4}{4!} + (-y_2) \frac{(-y_1)^2}{2!} +\frac{(-y_2)^2}{2!} + (-y_3)(-y_1) +(-y_4) \right)
$$
In this way you may generate the solution to any order. The expressions become too cumbersome to type here, but Mathematica has no trouble solving them. Here is a plot:

Since this is a perturbation series, the most we can hope for is that series is asymptotic. If you wished to take this further, note that after you have the $N$th order solution $y_N$, you have a sum of finitely many known terms
$$
y=\sum\limits_{n=0}^N \varepsilon^n y_n
$$
Then you could consider various methods of re-summing this series in $\varepsilon$. You could also investigate slightly different perturbation problems, eg: $y'''=e^{-\varepsilon y}$. I've tried the perturbation ansatz on the differential equation for $u(y)$, but that leads to more intractable differential equations at orders beyond first.
