#This worksheet illustrates using the Laplace transform to solve linear system
#of ODEs (and also just using Maple's dsolve command).

#Example System: Consider the linear constant-coefficient nonhomogeneous linear
#system of ODEs
var('t')
x1 = function('x1')(t);
x2 = function('x2')(t);
de1 = diff(x1,t) + 5*x1 - 6*x2 == -2*sin(t) + 10*cos(t) - 6*exp(t)
de2 = diff(x2,t) + 3*x1 - x2 == 6*cos(t)
#for functions x1(t) and x2(t), with initial data x1(0) = 2, x2(0) = 1.
#NOTE: The laplace command seems to be more successful when all dependent
#variables appear on the left side of the ODEs, all independent on the right.

#To solve, Laplace transform both sides of both equations
var('s');
de1lap = laplace(de1,t,s)
de2lap = laplace(de2,t,s)

#Substitute in the initial conditions, and (for convenience) let
#X1 = laplace(x1(t),t,s) and X2 = laplace(x2(t),t,s):
var('X1 X2');
de1lap2 = de1lap.substitute(x1(0)==2,x2(0)==1,laplace(x1(t),t,s)==X1,laplace(x2(t),t,s) == X2)
de2lap2 = de2lap.substitute(x1(0)==2,x2(0)==1,laplace(x1(t),t,s)==X1,laplace(x2(t),t,s) == X2)

#Solve for the transforms X1 and X2
sol = solve([de1lap2,de2lap2],[X1,X2])

#Assign transforms to variables X1sol and X2sol:
X1sol = X1.substitute(sol[0])
X2sol = X2.substitute(sol[0])

#Inverse transform to find the solutions
x1sol(t) = inverse_laplace(X1sol,s,t)
x2sol(t) = inverse_laplace(X2sol,s,t)
print("x1(t) = ",x1sol(t))
print("x2(t) = ",x2sol(t))

#Solution via dsolve: The solution can be obtained with Sage's dsolve_system()
#command:
sol = desolve_system([de1, de2], [x1,x2], ics=[0,2,1]);
print("Solution via desolve_system() is ",sol)