%% Variation of parameters for ODE system

clearvars

%% Example: \vec y\prime(t)= \left[\matrix{ 1 & 2 \cr 4 & 3 }\right]\vec y(t)+
% \left[\matrix{ e^{-t} \cr 1} \right]
%% Find fundamental set of solutions for homogeneous problem \vec y\prime=A\vec y
% Find eigenvalues and eigenvectors of matrix A. Find a fundamental set of solutions
%        \Psi(t) = \left[ \,\vec Y^{(1)}(t)\;,\vec Y^{(2)}(t) \,\right]
% See the class of April 10 for details. Here we use dsolve for the homogeneous problem:

A = sym([1 2; 4 3]);
syms y(t) [2 1]
syms C1 C2
sol = dsolve(diff(y)==A*y);
ysol = subs(y(t),sol)
Y1 = subs(ysol,{C1,C2},{1,0})
Y2 = subs(ysol,{C1,C2},{0,1})
Psi = [Y1,Y2]

%% Find \vec u\prime(t) by solving the linear system \Psi(t)\vec u\prime(t)=\vec f(t)

f = [exp(-t);1]

% Let \vec w(t) = \vec u\prime(t), this is the solution of the linear system \Psi(t) \vec w(t) = \vec f(t):

w = Psi\f          % solve linear system  Psi w = f 

% With the inverse matrix \Psi(t)^{-1}

inv(Psi)

% we can write the solution as w=\Psi(t)^{-1}\vec f(t)

w = inv(Psi)*f

%% Take the antiderivative to get \vec u(t)

u = int(w,t)

%% We obtain the particular solution \vec Y(t) = \Psi(t) \vec u(t)

Y = expand(Psi*u)

%% For comparison: Solve the inhomogeneous ODE with dsolve

sol = dsolve(diff(y)==A*y+f);
ysol = expand( subs(y(t),sol) )

% We get a particular solution by plugging in C_1=0,\;C_2=0:

Y = expand( subs(ysol,{C1,C2},{0,0}) )

% This is the same particular solution we obtained with variation of parameters above.