Solving ODEs with the Laplace Transform in Matlab

This approach works only for

linear differential equations with constant coefficients
right-hand side functions which are sums and products of
- polynomials
- exponential functions
- sine and cosine functions
- Heaviside (step) functions
- Dirac (impulse) ``functions''
initial conditions given at t = 0

The main advantage is that this approach allows to find symbolic solutions for right-hand side functions which are piecewise defined, and which contain Dirac impulse ``functions''.

You must first save the file Heaviside.m in your home directory. (This fixes a bug in Matlab.)

Example with piecewise defined right-hand side function

Consider the initial value problem

y'' + 3 y' + 2 y = f(t) , y(0) = 2 , y'(0) = 3

with the right-hand side function

f(t) = 1 for t<0
f(t) = t -2 for 3<t<6
f(t) = 2 for t>6

Define the necessary symbolic variables:

syms s t Y

As the right-hand side function is piecewise defined, rewrite it in terms of the Heaviside function:

f = '1 + (t-2-1)*Heaviside(t-3) + (2-(t-2))*Heaviside(t-6)'
ezplot(f,[0,10])

Find the Laplace transform of the right hand side function f(t):

F = laplace(f,t,s)

Find the Laplace transform of y'(t) : Y₁ = s Y - y(0)

Y1 = s*Y - 2

Find the Laplace transform of y''(t) : Y₂ = s Y₁ - y'(0)

Y2 = s*Y1 - 3

Set the Laplace transform of the left hand side minus the right hand side to zero and solve for Y:

Sol = solve(Y2+3*Y1+2*Y-F,Y)

Find the inverse Laplace transform of the solution:

sol = ilaplace(Sol,s,t)

Plot the solution:

ezplot(sol,[0,10])

Example with Dirac ``function''