#
# A Maple procedure is written to solve an ODE using variation of 
# parameters.  Thge proc. is tested on an ODE which yields a sol'n
# in terms of the Maple function Ei.  This solution is then analyzed
# graphically.
#
# procedure: variation_of_parameters
#  input:  two linearly indepent solutions to the homogeneous ODE
#          right hand side of ODE
#  output: solution to ODE
#
> variation_of_parameters:=proc(y1,y2,g) \
local y,u1,u2,y1p,y2p,w1,w2,wronski,x,t; \
y1p:=D(y1); \
y2p:=D(y2); \
wronski:=y1(t)*y2p(t)-y2(t)*y1p(t); \
wronski:=simplify(wronski); \
w1:=-y2(t)*g(t)/wronski; \
w2:=y1(t)*g(t)/wronski; \
u1:=int(w1,t=0..x); \
u2:=int(w2,t=0..x); \
y:=unapply(simplify(u1*y1(x)+u2*y2(x)),x); \
eval(y); \
end;
--------------------------------------------------------------------------------
> variation_of_parameters(cos,sin,sin);
--------------------------------------------------------------------------------
> 
# Now consider the problem y''+6*y'+8*y=1/(1+x^2). This is not at all 
# suited to Undetermined Coefficients. V of P is the only way to go. The two
# homogeneous solutions are y1:=exp(-4*x) and y2:=exp(-2*x). The V of P proc.
# we cooked up above takes functions as input and returns
# functions, not expressions, so we'll have to express y1 and y2 as functions:
> y1:=x->exp(-4*x);y2:=x->exp(-2*x);
--------------------------------------------------------------------------------
> 
# We also have to express g as a function:
> g:=x->1/(1+x^2);
# Now complicated DE's have complicated solutions, and no method will 
# ward off this reality. Let's see what happens. 
# 
--------------------------------------------------------------------------------
> yp:=variation_of_parameters(y1,y2,g);
# 
#
#  Maple Graphic.
#
# 
# Whoa. What ARE these functions? Well, ?Ei should tell us something.  
# Ei(n,x) = int(exp(-x*t)/t^n, t=1..infinity)   for  Re(x) > 0. 
# 
# In terms of material we shall be meeting shortly, Ei(n,x) is the 
# Laplace transform of the function u(t-1)/t^n. It is remarkable 
# how many things can be solved in closed form if one expands the list of 
# allowed functions to include the finite, though extensive, 
# list of special functions known to Maple. The point is, with a 
# closed-form solution there is a  good chance of making an analysis  
# of the long-range behavior of the solution. With numerical methods 
# alone, you're studying the solution with a flashlight, and you 
# never see anything but what you shine the beam of arithmetic on.
--------------------------------------------------------------------------------
> yp(0);
--------------------------------------------------------------------------------
> plot(yp(x),x=-1/4..2);
# 
--------------------------------------------------------------------------------
> 
# (You can also view the resulting plot directly from Mosaic.)
# 
# The solution tends slowly to zero as x increases. The reason is that 
# the solution would decay like e^(-2x) but for the constant  
# trickle of impetus it gets from the forcing function 1/(1+x^2). 
# Eventually the solution will settle down to close to (1/(8*x^2)). As
# x goes to -Infinity, the solution grows exponentially like exp(4*|x|). 
--------------------------------------------------------------------------------
> 
# Ei(n,x) = int(exp(-x*t)/t^n, t=1..infinity)   for  Re(x) > 0