One simple example of a nonlinear ODE with this property is given by

F(y, y', y", ...)=0, where F is a homogeneous function of its arguments, i.e.,

F(ax_1, ax_2, ...) = a^p F(x_1, x_2, ...)

for all a, x_1, x_2, etc. Then if y=y(x) is a solution, so is by(x) + cy(x), where b and c are arbitrary constants. The more interesting question is when the two solutions are *not* proportional and yet sum to give a new solution.

There is also a number of ODEs and systems of ODEs that admit *nonlinear* superposition of solutions. A certain broad class of such systems, known today as the Lie systems, was constructed already by Sophus Lie of the Lie group fame, see e.g. the papers by J.F. Carinena, J. de Lucas and M. Ranada http://iopscience.iop.org/1751-8121/41/30/304029/ and http://arxiv.org/abs/1103.4166 for a recent survey; see also the papers by R. L. Anderson (http://www.sciencedirect.com/science/article/pii/0167278982900586 and http://link.springer.com/chapter/10.1007%2F978-1-4684-4106-2_17) and by other authors (see e.g. T. C. Bountis, V. Papageorgiou and P. Winternitz http://scitation.aip.org/content/aip/journal/jmp/27/5/10.1063/1.527128).

Such (systems of) ODEs are intimately related to integrable (solitonic) nonlinear PDEs that admit nonlinear superposition of solutions, see e.g. this survey and references therein: P. Winternitz, Backlund transformations as nonlinear ordinary differential, or difference equations with superposition formulas. In: Backlund and Darboux transformations. The geometry of solitons (Halifax, NS, 1999), 429–436, CRM Proc. Lecture Notes, 29, Amer. Math. Soc., Providence, RI, 2001.

Well, there is a somewhat trivial answer to that question.

Suppose, your equation can be written in the general form

L[y] + N[y] = 0,

where L is some linear differential operator and N a non-linear differential operator. Then a sufficient condition for z = y_1+y_2 being a solution, provided y_1 and y_2 are solutions themselves is

N[y_1] + N[y_2] = N[y_1+y_2]. (A)

If we allow N to contain an additive free constant, then it is even sufficient, if we have

N[y_1] + N[y_2] = N[y_1+y_2] +c, (B)

where c is some new constant. This is obviously the case realized in the equation considered by you previously

y''' - y'^2 + x y - c = 0 (1)

where N[y] = - y'^2 - c.

One solution was y1 = -6/x, with c = -6, the other y2 = x^3/9, with c = 2/3. Then we have

N[y_1] + N[y_2] = -36/x^2 + 6 - x^4/9 - 2/3 = -36/x^2 - x^4/9 + 16/3 (*)

and

N[y_1+y_2] = -(6/x^2+x^2/3)^2 - c = -36/x^2 - 4 - x^4/9 - c (**)

(*) and (**) are identical if c = -28/3.

As you have observed yourself, this is due to the fact that the mixed term in N[y_1+y_2] becomes a constant.

I doubt that one can get more specific than (A) or (B) without additional restrictions on the form of the equations to be considered. [(A) holds for the original form of your equation, obtained by differentiating (1) once. Then N[y] = - 2 y' y''.]