A Substitution Based Method for Solving Non-Homogeneous Linear Differential Equations

1. Introduction

To solve higher-order non-homogeneous linear differential equations, the standard analytical frameworks are the Method of Undetermined Coefficients and the more general Variation of Parameters [1,2,3]. These approaches typically require first solving the associated homogeneous equation ($f\left(t\right)=0$) to obtain the complementary solution $y_h\left(t\right)$, which characterizes the system’s fundamental response. Subsequently, a particular solution $y_p\left(t\right)$ is sought to account for the non-homogeneous forcing function $f\left(t\right)$, such that the general solution is given by:

$$y\left(t\right)=y_h\left(t\right)+y_p\left(t\right)$$

In this work, an alternative approach based on operator factorization and sequential substitution is presented. By decomposing a higher-order linear operator into a product of first-order factors, the problem is reduced to a system of first-order linear equations. This method provides explicit integral representations for particular solutions without the need for heuristic guessing, Wronskian determinants, or transform methods.

2. The Homogeneous Case with Constant Coefficients

Consider the second-order homogeneous equation $y^{\prime \prime }+a{y}^{\prime }+by=0$ where $a,b\in \mathbb{R}$. In operator form, this is expressed as:

$$(D^2+aD+b)y=0$$

where $D={\displaystyle \frac{d}{dt}}$. The solution space is governed by the roots ${\lambda }_1,{\lambda }_2$ of the characteristic equation ${\lambda }^2+a\lambda +b=0$.

2.1. Real and Distinct Roots (${\lambda }_1\ne {\lambda }_2$)

Factoring the operator as $(D-{\lambda }_1)(D-{\lambda }_2)y=0$ and letting $u=(D-{\lambda }_2)y$, we obtain:

$$(D-{\lambda }_1)u=0\Rightarrow u\left(t\right)=c_1{e}^{{\lambda }_{1}t}$$

Substituting back into the definition of u:

$$y^{\prime }-{\lambda }_{2}y=c_1{e}^{{\lambda }_{1}t}$$

Applying the integrating factor $\mu \left(t\right)=e^{-{\lambda }_{2}t}$ yields:

$${\displaystyle \frac{d}{dt}}\left(e^{-{\lambda }_{2}t}y\right)=c_1{e}^{({\lambda }_1-{\lambda }_2)t}\Rightarrow y\left(t\right)=C_1{e}^{{\lambda }_{1}t}+C_2{e}^{{\lambda }_{2}t}$$

2.2. Repeated Real Roots (${\lambda }_1={\lambda }_2=\lambda$)

The operator ${(D-\lambda )}^{2}y=0$ with the substitution $u=(D-{\lambda }_2)y$ leads to $u^{\prime }-\lambda u=0$, giving $u=c_1{e}^{\lambda t}$. The subsequent equation $y^{\prime }-\lambda y=c_1{e}^{\lambda t}$ yields:

$${\displaystyle \frac{d}{dt}}\left(e^{-\lambda t}y\right)=C_1\Rightarrow y\left(t\right)=C_{1}t{e}^{\lambda t}+C_2{e}^{\lambda t}$$

2.3. Complex Conjugate Roots ($\lambda =a\pm bi$)

Factor as $(D-(a+bi\left)\right)(D-(a-bi\left)\right)y=0$ and let $u=(D-(a-bi\left)\right)y$. Then $u=C_1{e}^{(a+bi)t}$ where $C_1\in \mathbb{C}$. Solving for real y in the complex domain:

$$y={\displaystyle \frac{C_1}{2bi}}{e}^{(a+bi)t}+C_2{e}^{(a-bi)t}$$

By enforcing real-valued constraints $y=\overline{y}$ leads to $C_1=2bi{\overline{C}}_2$ and the standard form:

$$\begin{array}{cc}\hfill y\left(t\right)& =C_1{e}^{{\lambda }_{1}t}+{\overline{C}}_1{e}^{{\lambda }_{1}t}\hfill \\ \hfill & =e^{at}(c_{1}cos\left(bt\right)+c_{2}sin\left(bt\right))\hfill \end{array}$$

where $c_1,c_2\in \mathbb{R}$.

Consider $y^{\prime \prime }+a\left(t\right)y^{\prime }+b\left(t\right)y=0$. If the operator admits a factorization:

$$D^2+a\left(t\right)D+b\left(t\right)=(D-\alpha \left(t\right))(D-\beta \left(t\right))$$

where $a\left(t\right)=-(\alpha +\beta )$ and $b\left(t\right)=\alpha \beta -{\beta }^{\prime }$, then the solution is found sequentially. Let $u=(D-\beta (t\left)\right)y$; then:

$$u^{\prime }-\alpha \left(t\right)u=0\Rightarrow u\left(t\right)=C_1{e}^{\int \alpha \left(t\right)dt}$$

Therefore:

$$y^{\prime }-\beta \left(t\right)y=u\left(t\right)\Rightarrow y\left(t\right)=C_1{e}^{\int \beta \left(t\right)dt}\int e^{\int \left(\alpha \right(t)-\beta (t\left)\right)dt}dt+C_2{e}^{\int \beta \left(t\right)dt}$$

Example 2.1.

Consider the Cauchy–Euler equation

$$t^2{y}^{\prime \prime }-3t{y}^{\prime }+4y=0,t>0.$$

Observing that $t^2{D}^2-3tD+4=(tD-2)(tD-2)$, the equation becomes

$$(tD-2)(tD-2)y=0.$$

Let $u=(tD-2)y$, giving $t{u}^{\prime }-2u=0$ and $u=C_1{t}^2$ by separation of variables. Substituting $u=C_1{t}^2$ into $u=(tD-2)y$:

$$\begin{array}{cc}\hfill y^{\prime }-{\displaystyle \frac{2}{t}}y& =C_{1}t,\hfill \\ \hfill y& =C_1{t}^{2}ln\left|t\right|+C_2{t}^2.\hfill \end{array}$$

Remark. 

The standard form of $t^2{y}^{\prime \prime }-3t{y}^{\prime }+4y=0$ is $y^{\prime \prime }-{\displaystyle \frac{3}{t}}{y}^{\prime }+{\displaystyle \frac{4}{t^2}}y=0$, which can be represented as

$$\left(D^2-{\displaystyle \frac{3}{t}}D+{\displaystyle \frac{4}{t^2}}\right)y=\left(D-{\displaystyle \frac{1}{t}}\right)\left(D-{\displaystyle \frac{2}{t}}\right)y.$$

3. Non-Homogeneous Case with Constant Coefficients

For the non-homogeneous case $(D-{\lambda }_1)(D-{\lambda }_2)y=f\left(t\right)$, we apply the same substitution $u=(D-{\lambda }_2)y$. This yields the general integral representation of the particular solution $y_p$.

3.1. Real and Distinct Roots (${\lambda }_1\ne {\lambda }_2$)

Consider the non-homogeneous equation $(D-{\lambda }_1)(D-{\lambda }_2)y=f\left(t\right)$. By introducing the intermediate substitution $u=(D-{\lambda }_2)y$, we reduce the second-order problem to the following first-order system:

$$(D-{\lambda }_1)u=f\left(t\right)\Rightarrow u^{\prime }-{\lambda }_{1}u=f\left(t\right)$$

Solving via the integrating factor $\mu \left(t\right)=e^{-{\lambda }_{1}t}$, we obtain the general solution for $u\left(t\right)$:

$$u\left(t\right)=e^{{\lambda }_{1}t}\int e^{-{\lambda }_{1}t}f\left(t\right)dt+C_1{e}^{{\lambda }_{1}t}$$

Substituting this result back into the definition of u, we have:

$$y^{\prime }-{\lambda }_{2}y=e^{{\lambda }_{1}t}\int e^{-{\lambda }_{1}t}f\left(t\right)dt+c_1{e}^{{\lambda }_{1}t}$$

Applying the integrating factor $\mu \left(t\right)=e^{-{\lambda }_{2}t}$ to the second equation yields:

$$y\left(t\right)=e^{{\lambda }_{2}t}\int \left(e^{({\lambda }_1-{\lambda }_2)t}\int e^{-{\lambda }_{1}t}f\left(t\right)dt\right)dt+C_1{e}^{{\lambda }_{1}t}+C_2{e}^{{\lambda }_{2}t}$$

The particular solution $y_p\left(t\right)$ can be expressed in several equivalent forms. By symmetry of the operators $(D-{\lambda }_1)$ and $(D-{\lambda }_2)$, and by converting the nested indefinite integrals into a single definite integral over $[0,t]$, we derive the standard convolution form:

$$\begin{array}{cc}\hfill y_p\left(t\right)& =e^{{\lambda }_{2}t}\int \left(e^{({\lambda }_1-{\lambda }_2)t}\int e^{-{\lambda }_{1}t}f\left(t\right)dt\right)dt\hfill \\ \hfill & ={\int }_0^t{\displaystyle \frac{e^{{\lambda }_1(t-\tau )}-e^{{\lambda }_2(t-\tau )}}{{\lambda }_1-{\lambda }_2}}f\left(\tau \right)d\tau \hfill \end{array}$$

(1)

This representation explicitly reveals the Green’s function (kernel) for the operator:

$$G(t,\tau )={\displaystyle \frac{e^{{\lambda }_1(t-\tau )}-e^{{\lambda }_2(t-\tau )}}{{\lambda }_1-{\lambda }_2}}$$

3.2. Repeated Real Roots (${\lambda }_1={\lambda }_2=\lambda$)

In the case of a repeated characteristic root, the differential operator takes the form ${(D-\lambda )}^{2}y=f\left(t\right)$. By defining the intermediate variable $u=(D-\lambda )y$:

$$u^{\prime }-\lambda u=f\left(t\right)\Rightarrow u=e^{\lambda t}\int e^{-\lambda t}f\left(t\right)dt+C_1{e}^{\lambda t}$$

by the integrating factor method with $\mu \left(t\right)=e^{-\lambda t}$. Substituting this expression back into $u=(D-\lambda )y$, we have:

$$y^{\prime }-\lambda y=e^{\lambda t}\int e^{-\lambda t}f\left(t\right)dt+C_1{e}^{\lambda t}$$

Applying the same integrating factor $\mu \left(t\right)=e^{-\lambda t}$ to the second equation leads to:

$${\displaystyle \frac{d}{dt}}\left(e^{-\lambda t}y\right)=\int e^{-\lambda t}f\left(t\right)dt+C_1\Rightarrow y=e^{\lambda t}\int \left(\int e^{-\lambda t}f\left(t\right)dt\right)dt+C_{1}t{e}^{\lambda t}+C_2{e}^{\lambda t}$$

The particular solution $y_p\left(t\right)$ is characterized by the nested integral:

$$y_p\left(t\right)=e^{\lambda t}\int \left(\int e^{-\lambda t}f\left(t\right)dt\right)dt$$

Alternatively, using the identity for repeated integration, the particular solution can be expressed as a single definite integral, revealing the resonant kernel:

$$y_p={\int }_0^t(t-\tau )e^{\lambda (t-\tau )}f\left(\tau \right)d\tau$$

Example 3.1.

Consider the non-homogeneous equation $y^{\prime \prime }-2y^{\prime }+y=2e^t$. The characteristic equation ${\lambda }^2-2\lambda +1=0$ yields a repeated root $\lambda =1$. The associated differential operator is factored as ${(D-1)}^2$, and the forcing function is $f\left(t\right)=2e^t$. Utilizing the resonant kernel for repeated roots, the particular solution $y_p$ is given by:

$$y_p\left(t\right)=2{\int }_0^t(t-\tau )e^{(t-\tau )}{e}^{\tau }d\tau$$

Evaluating the integral with respect to τ yields:

$$y_p=t^2{e}^t.$$

3.3. Complex Conjugate Roots ($\lambda =a\pm bi$)

In the case of complex conjugate roots, the differential operator is factored as:

$$(D-(a+bi\left)\right)(D-(a-bi\left)\right)y=f\left(t\right).$$

Let ${\lambda }_1=a+bi$ and ${\lambda }_2=a-bi$. Then by Equation (1), the particular solution is given by:

$$y_p\left(t\right)={\int }_0^t{\displaystyle \frac{e^{{\lambda }_1(t-\tau )}-e^{{\lambda }_2(t-\tau )}}{{\lambda }_1-{\lambda }_2}}f\left(\tau \right)d\tau$$

To resolve this into a real-valued expression, we utilize Euler’s identity. The denominator simplifies to ${\lambda }_1-{\lambda }_2=2bi$, and the numerator becomes:

$$e^{(a+bi)(t-\tau )}-e^{(a-bi)(t-\tau )}=e^{a(t-\tau )}\left(e^{bi(t-\tau )}-e^{-bi(t-\tau )}\right)=2i{e}^{a(t-\tau )}sin\left(b(t-\tau )\right)$$

Substituting these into the integral, the $2i$ terms cancel, yielding the real-valued particular solution:

$$y_p\left(t\right)={\int }_0^t{e}^{a(t-\tau )}{\displaystyle \frac{sin\left(b\right(t-\tau \left)\right)}{b}}f\left(\tau \right)d\tau$$

Example 3.2.

Consider the non-homogeneous equation $y^{\prime \prime }-4y^{\prime }+5y=e^t$. The corresponding operator form is $(D-(2+i))(D-(2-i))y=e^t$, with $a=2$, $b=1$, and $f\left(t\right)=e^t$. Applying the formula derived above:

$$y_p\left(t\right)=e^{2t}{\int }_0^t{e}^{-2\tau }sin(t-\tau )e^{\tau }d\tau =e^{2t}{\int }_0^t{e}^{-\tau }sin(t-\tau )d\tau$$

Using the change of variables $u=t-\tau$ (where $du=-d\tau$), the integral transforms as follows:

$$y_p\left(t\right)=e^{2t}{\int }_t^0{e}^{-(t-u)}sin\left(u\right)(-du)=e^t{\int }_0^t{e}^{u}sin\left(u\right)du$$

Evaluation by integration by parts yields:

$$y_p\left(t\right)=e^t{\left[{\displaystyle \frac{e^u}{2}}(sinu-cosu)\right]}_0^t={\displaystyle \frac{e^{2t}}{2}}(sint-cost)+{\displaystyle \frac{1}{2}}{e}^t$$

Since the term ${\displaystyle \frac{e^{2t}}{2}}(sint-cost)$ is a linear combination of the complementary solutions, it may be absorbed into $y_h$. Thus, the simplest particular solution is:

$$y_p\left(t\right)={\displaystyle \frac{1}{2}}{e}^t$$

4. The Non-Homogeneous Case with Variable Coefficients

The substitution method generalizes to second-order linear differential equations with variable coefficients, provided the differential operator admits a factorization of the form:

$$L=D^2+a\left(t\right)D+b\left(t\right)=(D-\alpha \left(t\right))(D-\beta \left(t\right))$$

Expanding this operator product yields the following relations for the coefficient functions $\alpha \left(t\right)$ and $\beta \left(t\right)$:

$$\left\{\begin{array}{c}a\left(t\right)=-\left(\alpha \right(t)+\beta (t\left)\right)\hfill \\ b\left(t\right)=\alpha \left(t\right)\beta \left(t\right)-{\beta }^{\prime }\left(t\right)\hfill \end{array}\right.$$

To solve the non-homogeneous equation $(D-\alpha (t\left)\right)(D-\beta (t\left)\right)y=f\left(t\right)$, we define the intermediate variable $u=(D-\beta (t\left)\right)y$. This results in a first-order linear equation in u:

$$u^{\prime }-\alpha \left(t\right)u=f\left(t\right)$$

The solution for $u\left(t\right)$ is obtained via the standard integrating factor $e^{-\int \alpha \left(t\right)dt}$:

$$u\left(t\right)=e^{\int \alpha \left(t\right)dt}\int f\left(t\right)e^{-\int \alpha \left(t\right)dt}dt+C_1{e}^{\int \alpha \left(t\right)dt}$$

Substituting this expression back into the definition of u gives the first-order equation for y:

$$y^{\prime }-\beta \left(t\right)y=u\left(t\right)$$

Integrating once more, we arrive at the general solution $y\left(t\right)=y_p\left(t\right)+y_h\left(t\right)$, where the particular solution $y_p\left(t\right)$ is expressed by the nested integral formula:

$$y_p\left(t\right)=e^{\int \beta \left(t\right)dt}\int e^{\int \left(\alpha \right(t)-\beta (t\left)\right)dt}\left(\int f\left(t\right)e^{-\int \alpha \left(t\right)dt}dt\right)dt$$

Example 4.1.

Consider the non-homogeneous Cauchy–Euler equation $t^2{y}^{\prime \prime }-3t{y}^{\prime }+4y=2t^2$ for $t>0$. In standard form:

$$y^{\prime \prime }-{\displaystyle \frac{3}{t}}{y}^{\prime }+{\displaystyle \frac{4}{t^2}}y=2$$

We seek a factorization such that $D^2-{\displaystyle \frac{3}{t}}D+{\displaystyle \frac{4}{t^2}}=(D-\alpha \left(t\right))(D-\beta \left(t\right))$. Matching coefficients leads to the system:

$$\left\{\begin{array}{c}\alpha \left(t\right)+\beta \left(t\right)={\displaystyle \frac{3}{t}}\hfill \\ \alpha \left(t\right)\beta \left(t\right)-{\beta }^{\prime }\left(t\right)={\displaystyle \frac{4}{t^2}}\hfill \end{array}\right.$$

A solution to this system is $\alpha \left(t\right)={\displaystyle \frac{1}{t}}$ and $\beta \left(t\right)={\displaystyle \frac{2}{t}}$. Thus, the operator is factored as:

$$L=\left(D-{\displaystyle \frac{1}{t}}\right)\left(D-{\displaystyle \frac{2}{t}}\right)$$

The corresponding integral components are:

$$\left\{\begin{array}{cc}{e}^{-\int \alpha dt}\hfill & =e^{-lnt}=t^{-1}\hfill \\ e^{\int (\alpha -\beta )dt}\hfill & =e^{\int -t^{-1}dt}=t^{-1}\hfill \\ e^{\int \beta dt}\hfill & =e^{2lnt}=t^2\hfill \end{array}\right.$$

Substituting these components into the particular solution formula with $f\left(t\right)=2$:

$$y_p\left(t\right)=t^2\int {\displaystyle \frac{1}{t}}\left(\int {\displaystyle \frac{2}{t}}dt\right)dt$$

Successive integration yields:

$$y_p\left(t\right)=t^2\int {\displaystyle \frac{2lnt}{t}}dt=t^2{(lnt)}^2$$

5. Higher Order Differential Equations

The substitution method generalizes naturally to n-th order linear differential equations. If an n-th order linear differential operator L can be decomposed into a product of n first-order linear operators, the non-homogeneous problem is resolved by transforming it into a recursive chain of first-order equations. Consider the equation:

$$\left(D-{\alpha }_n\left(t\right)\right)\left(D-{\alpha }_{n-1}\left(t\right)\right)\dots \left(D-{\alpha }_1\left(t\right)\right)y=f\left(t\right)$$

where ${\alpha }_i\left(t\right)$ are real or complex-valued functions. We define a sequence of intermediate variables $\{u_1,u_2,\dots ,u_{n-1}\}$ as follows:

$$\begin{array}{cc}\hfill u_{n-1}& =(D-{\alpha }_{n-1})\dots (D-{\alpha }_1)y\hfill \\ \hfill u_{n-2}& =(D-{\alpha }_{n-2})\dots (D-{\alpha }_1)y\hfill \\ & ⋮\hfill \\ \hfill u_1& =(D-{\alpha }_1)y\hfill \end{array}$$

This substitution reduces the original n-th order equation into a system of n coupled first-order linear differential equations:

$$\begin{array}{cc}\hfill (D-{\alpha }_n)u_{n-1}& =f\left(t\right)\hfill \\ \hfill (D-{\alpha }_{n-1})u_{n-2}& =u_{n-1}\hfill \\ & ⋮\hfill \\ \hfill (D-{\alpha }_1)y& =u_1\hfill \end{array}$$

The general solution is obtained by solving each equation sequentially, starting from the top. Each step utilizes an integrating factor ${\mu }_i\left(t\right)=exp\left(-\int {\alpha }_i\left(t\right)dt\right)$. The particular solution $y_p\left(t\right)$ for the n-th order case can be expressed as a nested integral of depth n:

$$y_p\left(t\right)=e^{\int {\alpha }_{1}dt}\int e^{\int ({\alpha }_2-{\alpha }_1)dt}\int \dots \int e^{\int ({\alpha }_n-{\alpha }_{n-1})dt}\left(\int f\left(t\right)e^{-\int {\alpha }_{n}dt}dt\right){\left(dt\right)}^{n-1}$$

This hierarchical approach offers significant advantages over classical methods. Unlike the Method of Undetermined Coefficients, it requires no a priori assumption regarding the form of $f\left(t\right)$. Furthermore, it bypasses the computational burden of calculating the $n\times n$ Wronskian determinant required by the Method of Variation of Parameters, replacing it with n successive integrations.

6. Conclusions

The substitution-based method presented in this paper provides a robust and unified framework for the analysis of both homogeneous and non-homogeneous linear ordinary differential equations. By exploiting the algebraic structure of the differential operator, we transform a single n-th order problem into a recursive sequence of n first-order linear equations. For second-order ODEs, the particular solutions admit explicit integral representations as convolution kernels (Green’s functions), summarised in Table 1.

This systematic reduction offers several distinct advantages over classical pedagogical approaches. First, unlike the Method of Undetermined Coefficients, this framework requires no prior assumptions regarding the functional form of the forcing term $f\left(t\right)$. The particular solution $y_p$ emerges naturally through sequential integration, making the method applicable to a much broader class of non-elementary forcing functions that do not fit standard trial-solution templates.

Furthermore, the method bypasses the often prohibitive algebraic complexity inherent in the Method of Variation of Parameters. Specifically, it eliminates the necessity of constructing and evaluating $n\times n$ Wronskian determinants and performing subsequent matrix inversions. Instead, these are replaced by a more manageable chain of first-order integrating factors, which streamlines the calculation of the particular solution.

The approach also provides superior structural clarity in complex and resonant cases. The derivation of the sine-integral kernel for complex conjugate roots and the $(t-\tau )$ kernel for repeated roots arises directly from the operator factorization itself. This provides a clearer geometric and structural interpretation of resonance and oscillation than is typically found in standard transform-based or trial-solution methods. Moreover, as demonstrated with the Cauchy–Euler equation, the framework remains valid for variable-coefficient equations whenever a valid operator factorization $(D-\alpha (t\left)\right)(D-\beta (t\left)\right)$ can be identified, a process that often involves the resolution of an associated Riccati equation.

In summary, the substitution method bridges the gap between elementary solution techniques and more advanced operational calculus. By representing the particular solution as a nested integral—or a single integral via a reconstructed Green’s function—this approach offers both conceptual transparency and a powerful tool for exact symbolic integration in higher-order systems.