%% SECTION 4.3 - THE METHOD OF UNDETERMINED COEFFICIENTS

%%
% Here I solve the inhomogeneous fifth order ODE
%%
% 
% $$y^v  \,\mbox{--}\,3y^{iv} + 3y''' \,\mbox{--}\, 3y'' + 2y' = x^5 + 6e^{7x} 
% \,\mbox{--}\, 2e^{\,\mbox{--}\,x}\cos(3x) + 9 \sin(7x).$$
% 

%% Solving by hand (letting MATLAB do each calculation)
% First I have MATLAB carry out the same steps I would do if I were solving
% it by hand.  Then I see what happens if I ask MATLAB to solve the
% original equation.

%%
% Here is a MATLAB function (defined by a function-M file) which will apply
% the left side of the equation to a function f.

type diffop

%%
% The first step is to solve the corresponding homogeneous equation, by
% plugging 
% 
% $$y = e^{rx}$$
% 
% into the left side.

syms r x 
L = diffop(exp(r*x),x)


%% 
% I divide the left side by exp(rx).  I have to give the *simplify* 
% command to get MATLAB to carry out the division.

charpoly = simplify(L/exp(r*x))

%%
% Now I have to find the roots of the characteristic polynomial. *solve*
% assumes the right side of the equation is 0.

roots = solve(charpoly)

%%
% I know that 1, exp(x), exp(2x), cos(x) and sin(x) form a fundamental set
% of solutions of the homogeneous equation.

%%
% The right side of the equation has 4 terms, each of a different form, so 
% I break the inhomogeneous problem into 4 problems.  For the first, the
% right side is x^5.  So, since 1 solves the homogeneous equation, I look
% for a solution of the form y_1 equal to x times a polynomial of degree 5.  
% Here is a loop which builds it. 
y_1 = 0; 
for k = 0:5
vars{k+1} = ['a' num2str(k)];
u = vars{k+1}* x^(6-k);
y_1 = u + y_1;
end
y_1  

%%
% I plug it into the left side and subtract x^5.  The expression is very
% lengthy.  It helps to use *collect* to collect all the terms involving
% the same power of x (and factor the power out of the terms) and to use
% *pretty* so it doesn't run off the screen.

Y = collect(diffop(y_1,x) - x^5,x); pretty(Y)


%%
% I want to solve the system of equations that I get from setting the
% coefficents equal to 0.  I can either do this by copying and pasting the
% coefficients into the solve command or using a for loop to calculate the
% coefficients and set them equal to 0.  Here I use a loop to do it.  The
% value of the coefficient of x^j is the jth derivative of Y evaluated at
% 0. The corresponding equation is indexed by j+1.  
for k = 1:6 
    eqn = subs(diff(Y,x,k-1),x,0);
    eqns{k} = [char(eqn) '=0'];
end

%%
% Now I solve the resulting system.  
range = 1:6;
sol = solve(eqns{range},vars{range}); 
for k = range
    vals{k} = sol.(vars{k});
end


%%
% I substitute these into y_1.
Y_1 = subs(y_1,vars,vals)
%%
% I verify that this is a solution.  To get MATLAB to do the algebra after
% applying the left side to Y_1, I have to use *simplify*.

simplify(diffop(Y_1,x))


%%
% The right side of the second equation is 6 exp(7x), which doesn't solve 
% the homogeneous equation so I look for a solution of the form 
 
syms a

y_2 = a*exp(7*x)

%% 
% I plug it into the left side of the equation and subtract 6 exp(7x). 
% This time I want to collect the terms involving exp(7x).

Y = collect(diffop(y_2,x)-6*exp(7*x),exp(7*x))

%%
% I need the coefficient of exp(7x).

C = coeffs(Y,exp(7*x))

%%
% I solve this equation.  
[a] = solve(C)

%%
% I evaluate y_2 using this.
Y_2 = eval(y_2)

%%
% I verify that this is a solution.

diffop(Y_2,x)

%%
% The third term on the right is - 2exp(-x)cos(3x) which doesn't solve the 
% homoegeneous equation, so I look for a solution of the form 

clear a
syms a b
y_3 = a*exp(-x)*cos(3*x)+b*exp(-x)*sin(3*x)

%%
% I plug this into the left side of the equation and subtract - 2exp(-x)cos(3x)
% When I do this, every term will contain a factor of exp(-x), so 
% I divide the result by that.

Y = simplify((diffop(y_3,x) + 2*exp(-x)*cos(3*x))/exp(-x))

%%
% I could get the equations by evaluating Y and its derivative at 0.  
% Instead I will use  *coeffs* command.  I notice that Y is a linear
% "polynomial" in cos(3x) and sin(3x), so I want the coefficients of 
% cos(3x) and of sin(3x).

E1 = coeffs(Y,cos(3*x))
%%
% The coefficient of cos(3x) is the second term in E1
eq1 = E1(2)

%%
% Now I find the coefficient of sin(3x).

E2 = coeffs(Y,sin(3*x)); eq2 = E2(2) 

%%
% I solve the resulting equations using *solve*.  The *solve* command
% assumes that the right side of each equation is 0.  

[a,b] = solve(eq1,eq2); [a,b]

%%
% I evaluate y_3 using this.

Y_3 = eval(y_3)

%%
% I verify that this solves the equation.

diffop(Y_3,x)

%%
% The last term on the right is 9 sin(7x).  This doesn't solve the 
% homogeneous equation so I look for a solution of the 
% form

clear a b
syms a b
y_4 = a*cos(7*x) + b*sin(7*x)

%% 
% I plug this into the left side of the equation and subtract 9 sin(7x).

Y = diffop(y_4,x) - 9*sin(7*x)

%%
% Again I use *coeffs* to get the equations.

E1 = coeffs(Y,cos(7*x)); eq1 = E1(2)

%%
E2 = coeffs(Y,sin(7*x)); eq2 = E2(2)

%%
% I solve the equations.

[a,b] = solve(eq1,eq2)

%%
% The resulting solution of the equation is 

Y_4 = eval(y_4)

%%
% I verify that it solves the equation.

diffop(Y_4,x)
%%
% The resulting particular solution of the original equation is obtained by
% adding these 4 solutions.

y_p = Y_1 + Y_2 + Y_3 + Y_4; pretty(y_p)

%% Solving with Matlab

%%
% What happens if you let MATLAB solve the equation?
eqn =  'D5y - 3*D4y + 3*D3y - 3*D2y + 2*Dy= x^5 + 6*exp(7*x) - 2*exp(-x)*cos(3*x) + 9*sin(7*x)';
Y = dsolve(eqn, 'x');
pretty(simple(Y))
   
%%
% This is a general solution of the equation.  I compare it with our
% particular solution.

simplify(Y - y_p)

%%
% I get a solution of the homogeneous equation, which is what I should
% get.