Closed-form solution for a second order ODE I have a series of second order ODE's that are in the form:
$$
Ay''+By'+Cy=\frac{1}{1+e^{-t}}
$$
($y'' = \frac{d^2}{dt^2}y(t)$,  $y' = \frac{d}{dt}y(t)$ )
The initial conditions are known: $y(0)=y_0$, $y'(0)=y'_0$ 
Clearly, the general solution for the homogeneous form $Ay''+By'+Cy=0$ is trivial; I'm stuck in finding the particular solution.
Is there a hope that I can come up with a closed-form solution parametrized by $A$, $B$, $C$ and the initial conditions (All in $\mathbb{R}$)?
Background: Initially I was trying to solve the ODE with the step function (the heavy-side function: $u(t)=1~if~t>0; ~0~ow$) on the right-hand side and a different set of coefficients on the left-hand side:
$$
(A_1t+A2)y''+(B_1+B2)y'+(C_1t+C2)y=\frac{1}{1+e^{-t}}
$$
But, I thought maybe I could find the solution with a smooth function on the RHS, and a simpler ODE in the first place.
Therefore, a solution to the original (either with affine or with constant coefficients) problem will be helpful too.
Thanks for any help!
 A: As hinted, variation of parameters is one way to go. The homogeneous solution is given by
$$y_h(t)=\underbrace{c_1 e^{\lambda_1 t}}_{=:y_1(t)}+\underbrace{c_2 e^{\lambda_2 t}}_{=:y_2(t)}$$
where 
$$\lambda_1=\frac{-B-\sqrt{B^2-4AC}}{2A}, ~\lambda_2=\frac{-B+\sqrt{B^2-4AC}}{2A}$$
Now compute the Wronskian
$$\mathcal{W}(t)=\begin{vmatrix} y_1 & y_2 \\ \tfrac{d}{dt}y_1 & \tfrac{d}{dt}y_2  \end{vmatrix}=...=\frac{e^{-\frac{Bt}{A}}\sqrt{B^2-4AC}}{A} $$
Now the particular solution is given by
$$y_p(t)=-y_1(t) \int \frac{y_2(t)}{A(e^{-t}+1)\mathcal{W}(t)} ~\text{d}t +y_2(t) \int \frac{y_1(t)}{A(e^{-t}+1)\mathcal{W}(t)} ~\text{d}t$$
and the solution to the ODE is
$$y(t)=y_h(t)+y_p(t).$$
You asked for a closed form solution but this is quite strenuous since it involves hypergeometric functions. WolframAlpha and Mathematica tell me the following

UPDATE: I checked the system
$$(a_1t+a_2)y''(t)+(b_1t+b_2)y'(t)+(c_1t+c_2)y(t)=0$$
and according to WolframAlpha the solution is everything but nice; it involves the Laguerre polynomial and the confluent hypergeometric function of second kind. See yourself:

Maybe some numerical methods would be more appropriate. 
A: $$ay''+by'+cy=\frac{1}{1+{{e}^{-t}}}$$
Take the Laplace transform valid for t>0, 
$$a\left( -y{{'}_{0}}-s{{y}_{0}}+{{s}^{2}}Y \right)+b\left( sY-{{y}_{0}} \right)+cY=\int\limits_{0}^{\infty }{\frac{{{e}^{-st}}}{1+{{e}^{-t}}}dt}$$
so
$$Y\left( s \right)=\frac{a{{y}_{0}}s+\left( {{y}_{0}}b+ay{{'}_{0}} \right)}{a{{s}^{2}}+bs+c}+\frac{1}{a{{s}^{2}}+bs+c}\int\limits_{0}^{\infty }{\frac{{{e}^{-st}}}{1+{{e}^{-t}}}dt}$$
Note:
$$\int\limits_{0}^{\infty }{\frac{{{e}^{-st}}}{1+{{e}^{-t}}}dt}=\sum\limits_{n=0}^{\infty }{{{\left( -1 \right)}^{n}}}\int\limits_{0}^{\infty }{{{e}^{-\left( s+n \right)t}}dt}=\sum\limits_{n=0}^{\infty }{\frac{{{\left( -1 \right)}^{n}}}{s+n}}$$
And therefore upon inverting we have 
$$y\left( t \right)=\frac{1}{2\pi i}\int\limits_{{}}^{{}}{{}}\frac{a{{y}_{0}}s+\left( {{y}_{0}}b+ay{{'}_{0}} \right)}{a{{s}^{2}}+bs+c}{{e}^{st}}+\frac{1}{a{{s}^{2}}+bs+c}\sum\limits_{n=0}^{\infty }{\frac{{{\left( -1 \right)}^{n}}}{s+n}{{e}^{st}}}$$
Where the integral is the Bromwich contour. Let $a{{s}^{2}}+bs+c=a\left( s-{{\beta }_{+}} \right)\left( s-{{\beta }_{-}} \right)$ and let ${{\beta }_{\pm }}=\frac{-b\pm \sqrt{\Delta }}{2a}\Rightarrow a\left( {{\beta }_{+}}-{{\beta }_{-}} \right)=\sqrt{\Delta }\\$ where we have the discriminant $\Delta ={{b}^{2}}-4ac$.  Then
$$y\left( t \right)=\frac{1}{2\pi i}\int\limits_{{}}^{{}}{{}}\frac{a{{y}_{0}}s+\left( {{y}_{0}}b+ay{{'}_{0}} \right)}{a\left( s-{{\beta }_{+}} \right)\left( s-{{\beta }_{-}} \right)}{{e}^{st}}+\frac{1}{a\left( s-{{\beta }_{+}} \right)\left( s-{{\beta }_{-}} \right)}\sum\limits_{n=0}^{\infty }{\frac{{{\left( -1 \right)}^{n}}}{s+n}{{e}^{st}}}$$
Calculating residues we have
$$y\left( t \right)=\frac{1}{\sqrt{\Delta }}\left( a{{y}_{0}}{{\beta }_{+}}+\left( {{y}_{0}}b+ay{{'}_{0}} \right)+\sum\limits_{n=0}^{\infty }{\frac{{{\left( -1 \right)}^{n}}}{{{\beta }_{+}}+n}} \right){{e}^{{{\beta }_{+}}t}}-\frac{1}{\sqrt{\Delta }}\left( a{{y}_{0}}{{\beta }_{-}}+\left( {{y}_{0}}b+ay{{'}_{0}} \right)+\sum\limits_{n=0}^{\infty }{\frac{{{\left( -1 \right)}^{n}}}{{{\beta }_{-}}+n}} \right){{e}^{{{\beta }_{-}}t}}+\sum\limits_{n=0}^{\infty }{\frac{{{\left( -1 \right)}^{n}}}{a\left( n+{{\beta }_{+}} \right)\left( n+{{\beta }_{-}} \right)}{{e}^{-nt}}}$$
It is then a reasonably simple matter to show that this solution satisfies the DE and the boundary conditions.  For example, consider 
$$ay''+by'+cy=\frac{1}{\sqrt{\Delta }}\left\{ \begin{align}
  & \left( a{{\beta }^{2}}_{+}+b{{\beta }_{+}}+c \right)\left( a{{y}_{0}}{{\beta }_{+}}+\left( {{y}_{0}}b+ay{{'}_{0}} \right)+\sum\limits_{n=0}^{\infty }{\frac{{{\left( -1 \right)}^{n}}}{{{\beta }_{+}}+n}} \right){{e}^{{{\beta }_{+}}t}} \\ 
 & -\left( a{{\beta }^{2}}_{-}+b{{\beta }_{-}}+c \right)\left( a{{y}_{0}}{{\beta }_{-}}+\left( {{y}_{0}}b+ay{{'}_{0}} \right)+\sum\limits_{n=0}^{\infty }{\frac{{{\left( -1 \right)}^{n}}}{{{\beta }_{-}}+n}} \right){{e}^{{{\beta }_{-}}t}} \\ 
\end{align} \right\}+\sum\limits_{n=0}^{\infty }{\frac{{{\left( -1 \right)}^{n}}\left( a{{n}^{2}}-bn+c \right)}{a\left( n+{{\beta }_{+}} \right)\left( n+{{\beta }_{-}} \right)}{{e}^{-nt}}}$$
Now the first two terms disappear because $\beta $  is a solution of the quadratic.  We have then
$$ay''+by'+cy=\sum\limits_{n=0}^{\infty }{\frac{{{\left( -1 \right)}^{n}}\left( a{{n}^{2}}-bn+c \right)}{a\left( n+{{\beta }_{+}} \right)\left( n+{{\beta }_{-}} \right)}{{e}^{-nt}}}$$
Note  $a{{s}^{2}}+bs+c=a\left( s-{{\beta }_{+}} \right)\left( s-{{\beta }_{-}} \right)$ hence $a{{n}^{2}}-bn+c=a\left( -n-{{\beta }_{+}} \right)\left( -n-{{\beta }_{-}} \right)=a\left( n+{{\beta }_{+}} \right)\left( n+{{\beta }_{-}} \right)$. 
$$ay''+by'+cy=\sum\limits_{n=0}^{\infty }{{{\left( -1 \right)}^{n}}{{e}^{-nt}}}=\frac{1}{1+{{e}^{-t}}}$$
To check the boundary conditions again observe $a\left( {{\beta }_{+}}-{{\beta }_{-}} \right)=\sqrt{\Delta }$, but also $\left( {{\beta }_{+}}+{{\beta }_{-}} \right)=-b/a$.   Therefore after some simplification
$$y\left( 0 \right)=\frac{1}{\sqrt{\Delta }}\left\{ a{{y}_{0}}\left( {{\beta }_{+}}-{{\beta }_{-}} \right)+\sum\limits_{n=0}^{\infty }{\frac{{{\left( -1 \right)}^{n}}\left( {{\beta }_{-}}-{{\beta }_{+}} \right)}{\left( {{\beta }_{+}}+n \right)\left( {{\beta }_{-}}+n \right)}}+\sqrt{\Delta }\sum\limits_{n=0}^{\infty }{\frac{{{\left( -1 \right)}^{n}}}{a\left( n+{{\beta }_{+}} \right)\left( n+{{\beta }_{-}} \right)}} \right\}={{y}_{0}}$$
$$y'\left( 0 \right)=\frac{1}{\sqrt{\Delta }}\left\{ -\frac{b\sqrt{\Delta }}{a}{{y}_{0}}+\frac{\sqrt{\Delta }}{a}\left( {{y}_{0}}b+ay{{'}_{0}} \right)+\sum\limits_{n=0}^{\infty }{{{\left( -1 \right)}^{n}}\frac{\left( {{\beta }_{+}}-{{\beta }_{-}} \right)n}{\left( {{\beta }_{+}}+n \right)\left( {{\beta }_{-}}+n \right)}}-\sqrt{\Delta }\sum\limits_{n=0}^{\infty }{\frac{{{\left( -1 \right)}^{n}}n}{a\left( n+{{\beta }_{+}} \right)\left( n+{{\beta }_{-}} \right)}} \right\}=y{{'}_{0}}$$.
Now consider once more the solution y(t), and note it may be written in terms of some special functions
$$\begin{align}
  & y\left( t \right)=\frac{1}{\sqrt{\Delta }}\left( a{{y}_{0}}{{\beta }_{+}}+\left( {{y}_{0}}b+ay{{'}_{0}} \right)+\Phi \left( -1,1,{{\beta }_{+}} \right) \right){{e}^{{{\beta }_{+}}t}}-\frac{1}{\sqrt{\Delta }}\left( a{{y}_{0}}{{\beta }_{-}}+\left( {{y}_{0}}b+ay{{'}_{0}} \right)+\Phi \left( -1,1,{{\beta }_{-}} \right) \right){{e}^{{{\beta }_{-}}t}} \\ 
 & +\frac{1}{{{\beta }_{+}}{{\beta }_{-}}\sqrt{\Delta }}\left( {{\beta }_{+}}{}_{2}{{F}_{1}}\left( 1,{{\beta }_{-}};1+{{\beta }_{-}},{{e}^{-t}} \right)-{{\beta }_{-}}{}_{2}{{F}_{1}}\left( 1,{{\beta }_{+}};1+{{\beta }_{+}},{{e}^{-t}} \right) \right) \\ 
\end{align}$$
Where we have the Lerch transcendent and the hypergeometric function.  This may be of interest perhaps, if you are considering certain types of zeros of the quadratic.  For example note that $\Phi \left( -1,1,\beta  \right)$ becomes undefined whenever $\beta <0$ and an integer, which implies if we are considering two real integer solutions, then they must be positive which then places restrictions upon a,b,c etc.  Perhaps also of interest is that it takes on specific forms depending upon the roots of the quadratic (as does the hypergeometric function).  
