Procedure for solving higher order constant coefficient linear differential equations

Part 0. The equation:

Suibstitute your data for the following:

order  =  4;
coefficients  =  {1,  3,  4, 2,  
          5};
f[t_] := (E^t) (t^2)

It is required to solve the equation

\(eqn\  = \ \((\((∑\+\(j = 0\)\%order 
                  coefficients[[
                  j + 1]]\ ∂\_{t, j}y[
                    t])\) == f[t])\)\)

Part 1. Solve the characteristic equation

Consider the chartacteristic equation:

charEqn = (∑_(j = 0)^order) coefficients[[
              j + 1]] (r^j) == 0

charSoln  =  NSolve[charEqn,  r]

α  =  r  /.  charSoln

Part 2.  Solve the homogeneous equation

homogEqn  =  (a y''[t]  +  b y'[t]  +  
            c y[t]  ==  0)

A fundamental family of solutions to the homogeneous equation  depends on the nature of the roots of the characteristic equation.  The roots may appear in complex conjugate pairs and may have various multiplicities.

A root α with multiplicity one corresponds to a single solution (e^(α t)).

A root α with higher multiplicity  m corresponds to solutions 
(e^(α t)), t (e^(α t)), (t^2) (e^(α t)), ... (t^(m - 1)) (e^(α t)).

A complex conjugate pair  λ  + i μ,  λ  - i μ  corresponds to solutions (e^(λ t)) Cos[μ t], (e^(λ t))Sin[μ t].

In  our example, we have two complex conjugate pairs.

A fundamental family of solutions to the homogenous equation is:

φ = {(E^(Re[α[[1]]] t)) Cos[
            Im[α[[1]]] t], (E^(Re[α[[1]]] t)) Sin[Im[α[[1]]] t], (E^(Re[α[[3]]] t)) Cos[Im[α[[3]]] t], (E^(Re[α[[3]]] t)) Sin[Im[α[[3]]] t]}

The general solution to the homogenous equation is

(∑_(i = 1)^order) C[i] φ[[
          i]]

Part 2.  Find a particular solution to the inhomogenous equation.

We guess a  particular solution to the inhomogeneous equation of the following form.  (Substitute the appropriate functional form depending on the nature of the right hand side f[t].)

Clear[yp, A0,  A1,  A2]
yp[t_] :=  (E^t) (A0 + A1 t + A2 (t^2))

We put this proposed solution into the left hand side of the differential equation:

\(LHS\  = \ 
      Expand[\((∑\+\(j = 0\)\%order coefficients[[
                j + 1]]\ ∂\_{t, j}yp[
                  t])\)]\)

We set the left hand side equal to the right hand side, namely f[t]

coeffEqn  =  (LHS  ==  f[t])

Substituting n different values of t yields n linear equations in the n unknown parameters  A0, A1, A2, ....

linearEquations  =  {
(coeffEqn1  =  
            coeffEqn  /.  
              t ->  0),  
(coeffEqn2  =  
            coeffEqn  /.  
              t   ->  1),  
(coeffEqn3  =  
            coeffEqn  /.  t  ->  2)}

coeffSoln  =  Solve[linearEquations,  {A0,  A1,  A2}]

YP[t_]  :=  (yp[t] /.  
          coeffSoln[[1]])
YP[t]

Part 3.  The general solution

The general solution to the differential equation is:

genSoln[
        t_]  :=  (∑_(i = 1)^order) CC[
              i] φ[[i]]  +  
        YP[t]
genSoln[t]

Part 4.  An initial value problem.

Let's consider the IVP consisting of our differential equation

Clear[y]
\(eqn\  = \ \((\((∑\+\(j = 0\)\%order 
                  coefficients[[
                  j + 1]]\ ∂\_{t, j}y[
                    t])\) == f[t])\)\)

with the initial conditions:

Clear[y, InitialConditions]
InitialConditions={y[0] == 3,y'[0]==2,y''[0]==0,y'''[0] == 4}

If we substitute the general solution for y[t], the initial conditions become 2 linear equations in the unknowns C1 and C2.

y[t_]  =  
        genSoln[t];

    InitialConditions

We solve these linear equations and substitute the result into y[t] to obtain the solution to the IVP.

IVSoln  =  Solve[
          InitialConditions,  {CC[1],  CC[2],  CC[3],  
            CC[4]}][[1]]
F[t_]  =  
        Expand[(y[t]  /.  IVSoln)];
F[t]

Plot[F[t],  {t,  -5,  5}]

Part 5. Compare with built in Mathematica routines

Clear[y,  F]
\(eqn\  = \ \((\((∑\+\(j = 0\)\%order 
                  coefficients[[
                  j + 1]]\ ∂\_{t, j}y[
                    t])\) == f[t])\)\)
dSoln  =  DSolve[eqn,  y[t],  t];
F[t_]  =  Expand[ Simplify[y[t] /.  dSoln[[1]]]]