Inhomogeneous 2nd-order linear differential equation I need to solve this:
$y'' + ay' + by = x(t)$
where nothing about the form of $x$ is known, except that it is bounded and non-negative.  In addition it is known that $y(0) = 0$ and $y'(0) = 0$ (and $y''(0) = 0$ as well). I plugged this into the Wolfram ODE solver and got back a 'possible Lagrangian', which I really don't know how to use.  (I have an undergraduate degree in math but I never took differential equations.)  I also found this resource pt1, pt2, but it assumes that the solution is exponential in nature, which seems to violate my initial conditions since exponentials are never 0.  I'm hoping there's a way to get my head around this without putting myself through the full course of diff eq so any help would be much appreciated.
Update: So actually Wolfram Alpha will solve this equation (simplified slightly in this edit), it was just the 'widget' that stopped short.  However, it, as well as the techniques given in the references, explicitly assumes that the solution has an exponential form.  Apparently this is a standard assumption but it is not the case for my data.  $y$ is roughly sigmoidal.  Coding up Wolfram's solution and choosing arbitrary constants $c_1$, $c_2$ so as to minimize the error in a case where $y$ was known produced a clearly wrong answer.
So, I will sharpen my question: How to solve the above differential question without assuming that $y$ is proportional to $e^{\lambda t}$?
 A: Let $f_\pm (t) := \exp(\lambda_\pm t)$ where $\lambda_\pm := (-a \pm \sqrt{a^2 - 4 b})/2$ be the fundamental solutions to the homogenous equation. Then the solution to the ODE, subject to $y(0)=r_0$ and $y^{'}(0)= \mu_0$, reads:
\begin{equation}
y(t) = \left(\begin{array}{cc} f_+(t) & f_-(t) \end{array}\right) \cdot
\left( \begin{array}{cc} f_+(0) & f_-(0) \\ f_+^{'}(0) & f_-^{'}(0) \end{array} \right)^{-1} \cdot
\left( \begin{array}{c} r_0 \\ \mu_0 \end{array} \right)
+
\int\limits_0^t \left(\frac{f_+(\xi) f_-(t)-f_-(\xi) f_+(t)}{f_+(\xi) f_-^{'}(\xi) - f_-(\xi) f_+^{'}(\xi)}\right) x(\xi) d\xi
\end{equation}
This can be proven using the Green's function method, for instance.
 Now, the first term on the right hand side is the generic solution to the homogenous equation and the second term is a specific solution to the inhomogenous equation.The quantity in the denominator in the integrand is called the Wronskian $W(\xi)$ and it reads:
\begin{equation}
W(\xi) = (\lambda_- - \lambda_+) e^{(\lambda_- + \lambda_+) \xi} = - \sqrt{a^2 - 4 b} e^{-a \xi}
\end{equation}
 In our case $r_0=0$ and $\mu_0 = 0$ and so the solution reads:
\begin{equation}
y(t) = \frac{-1}{\sqrt{a^2 - 4 b}} 
\left( 
e^{\lambda_- t} \int\limits_0^t e^{(\lambda_++a) \xi} x(\xi) d \xi
-
e^{\lambda_+ t} \int\limits_0^t e^{(\lambda_-+a) \xi} x(\xi) d \xi
\right)
\end{equation}
A: It's simple with Operator Calculus .
Quite analogous to Example 1 on page 3 of this PDF we have:
$$ y'' + a y' + b y = x(t)$$
Operator Calculus enables us to "abstract" as much as possible from the solution $y(t)$ :
$$ \left[ \left(\frac{d}{dt}\right)^2 + a \frac{d}{dt} + b \right] y(t) = x(t)$$
What we are going to do next is decompose into 
factors. We are (almost) forced to do so, because the operator $(d/dt)$ and
the constants $a$ and $b$ mutually behave as if they were ordinary
numbers: the commutator of a differentiation and a constant is zero. But
the commutative law for ( linear ) operators is the only thing which could
be deviant from algebra with ordinary numbers.
$$ \left(\frac{d}{dt} \right)^2 + a \left(\frac{d}{dt} \right) + b =
   \left(\frac{d}{dt} - \lambda_1 \right)
   \left(\frac{d}{dt} - \lambda_2 \right) $$
$\lambda_1$ and $\lambda_2$ are roots of the so-called characteristic
equation : $$ \lambda^2 + a \lambda + b = 0 \quad \Longrightarrow \quad \lambda_{12}= (-a \pm \sqrt{a^2 - 4 b})/2 $$
Herewith, the differential equation can be rewritten as:
$$ \left(\frac{d}{dt} - \lambda_1 \right)
   \left(\frac{d}{dt} - \lambda_2 \right) y(t) = x(t) $$
We are going to use now the "very useful formula" from Operator Calculus
(please find it on page 2 of the abovementioned document) :
$$
   \frac{d}{dt} + f = e^{-\int f\,dt}\; \frac{d}{dt} \; e^{+\int f\,dt}
$$
In our case:
$$ \frac{d}{dt} - \lambda_{1,2} =
   e^{ - \int - \lambda_{1,2} \, dt} \ \frac{d}{dt}\
   e^{ + \int - \lambda_{1,2} \, dt} =
 e^{\, \lambda_{1,2} t } \ \frac{d}{dt}\ e^{\, - \lambda_{1,2} t } $$
Giving, at last, for the O.D.E. :
$$ e^{ \lambda_1 t } \frac{d}{dt}\ e^{ - \lambda_1 t }\;
   e^{ \lambda_2 t } \frac{d}{dt}\ e^{ - \lambda_2 t } \; y(t) = x(t) $$
Systematic integration is possible now:
$$ e^{-\lambda_1 t} e^{\lambda_2 t} \frac{d}{dt}\ e^{-\lambda_2 t} y(t) =
  \int x(t) e^{ - \lambda_1 t }\,dt $$
$$ y(t) = e^{\lambda_2 t} \int e^{(\lambda_1-\lambda_2) t} 
   \left[ \int x(t) e^{ - \lambda_1 t }\,dt \right] \, dt $$
As has been said, $\lambda$ can be solved from $\lambda^2 + a\lambda + b = 0$,
a quadratic equation with discriminant: $ D=a^2-4b$.
Real valued as well as complex solutions are possible. A special case for $D=0$
may be distinguished; in this case $ \lambda_1 = \lambda_2 $
and $e^{(\lambda_1-\lambda_2)t} = 1$ .
EDIT. Almost forgot the boundary conditions $y(0)=y'(0)=0$ to take into account:
$$ y(t) = e^{\lambda_2 t} \int_0^t e^{(\lambda_1-\lambda_2) u} 
   \left[ \int_0^u x(v) e^{ - \lambda_1 v }\,dv \right] \, du $$
Note. If "nothing about the form of x is known" - but why a bounded & non-negative requirement - then $y''(0) = 0$ is in general not true , unless $x(0) = 0$ :
$$
y''(t) = \lambda_2^2 e^{\lambda_2 t} \int_0^t e^{(\lambda_1-\lambda_2) u} 
   \left[ \int_0^u x(v) e^{ - \lambda_1 v }\,dv \right] \, du +
   (\lambda_1 + \lambda_2) e^{\lambda_1 t} \int_0^t x(v) e^{ - \lambda_1 v }\,dv
   \; + \; x(t)
$$
With differential equations of the kind, $\,y''(0) = 0\,$ cannot serve as a boundary condition anyway.
A: I use the Laplace transform (see e.g. here ).
Taking the Laplace transform of both sides we obtain
$$
(as^2+bs+c)Y(s)=X(s),
$$
$Y(s),X(s)$ are the Laplace transform of $y(t),x(t)$, respectively. (In fact, $x(t)$ can be distribution also, the so-called generalized function.)
So we have 
$$
Y(s)=\frac{X(s)}{as^2+bs+c}.
$$
Using the convolution theorem we obtain 
$$
y(t)=x(t)*\mathcal{L}^{-1}\left(\frac{1}{as^2+bs+c} \right),
$$
where $(f*g)(t)=\int_0^t f(t-v)g(v)\,dv$, and $f,g$ are piecewise continuous on $[0,\infty)$ (but this assumption can be relaxed). 
Using the partial fraction decomposition 
$$
\frac{1}{as^2+bs+c}=\frac{1}{a}\left(\frac{c_1}{s-s_1}+\frac{c_2}{s-s_2} \right),
$$
or
$$
\frac{1}{as^2+bs+c}=\frac{1}{a(s-s_1)^2}.
$$
Here
$$
\mathcal{L}^{-1}\left(\frac{1}{s-A} \right)=e^{At}, \qquad\mathcal{L}^{-1}\left(\frac{1}{(s-A)^2} \right)=e^{At}t.
$$
A: How about using variation of parameters? 
If you name $y_i(t)$ as the two parts of the homogeneous solution, obtained by assuming exponentials, so $$y(t) = A y_1(t) + B y_2(t) + y_p(t),$$ where $y_p(t)$ is the particular solution of the ODE, then by using variation of parameters as follows:
$$y(t) = A(t) y_1(t),$$
you will conclude that:
$$y_p(t) = y_1(t) \int M(t) \, dt, $$
where $M(t)$ is given by the following expression:
$$M(t) = \frac{e^{-at}}{y_1^2(t)} \int e^{at} y_1(t) x(t) \, dt.$$
You should apply now the initial conditions to the final solution, $y(t)$. 
If you want me to elaborate some more on my answer, please let me know.
Cheers!
