#
# Hand technique of variation of parameters is illustrated using Maple.
#
# VARIATION of PARAMETERS
# Dr. Kirschner, Math 308.
# 
> restart;
> diffeq:=diff(y(x),x$2)+2*diff(y(x),x)+y(x)=3*exp(-x);

                     /   2      \
                     |  d       |     /  d      \
           diffeq := |----- y(x)| + 2 |---- y(x)| + y(x) = 3 exp(- x)
                     |   2      |     \ dx      /
                     \ dx       /
> char_eq:=r^2+2*r+1; solve(",r);

                                         2
                             char_eq := r  + 2 r + 1

                                     -1, -1
> with(linalg): 
Warning: new definition for   norm
Warning: new definition for   trace

> FS:=[exp(-x),x*exp(-x)]; WW:=Wronskian(FS,x);

                          FS := [exp(- x), x exp(- x)]

                         [  exp(- x)         x exp(- x)      ]
                   WW := [                                   ]
                         [ - exp(- x)  exp(- x) - x exp(- x) ]
> R:=matrix(2,1,[0,3*exp(-x)]);\


                                    [      0     ]
                               R := [            ]
                                    [ 3 exp(- x) ]
# MATRIX MULTIPLICATION IS WITH A &
# 
> (1/WW) &* R; evalm(");

                                      1
                                    ---- &* R
                                     WW

                                    [ - 3 x ]
                                    [       ]
                                    [   3   ]
# NOW< INTEGRATE BOTH SIDES TO GET u,v
# (map allows us to do this to both sides at once)
#  WE will use the constants so when we are done
# we will get both the homogeneous solution and the
# particular solution at the same time!
> map(int,",x)+matrix(2,1,[C1,C2]); U:=evalm(");

                              [        2 ]
                              [ - 3/2 x  ]   [ C1 ]
                              [          ] + [    ]
                              [    3 x   ]   [ C2 ]

                                  [        2      ]
                                  [ - 3/2 x  + C1 ]
                             U := [               ]
                                  [    3 x + C2   ]
> sol:=y(x)=dotprod(FS, convert(U,vector)); 

                                        2
         sol := y(x) = exp(- x) (- 3/2 x  + C1) + x exp(- x) (3 x + C2)
> C1:=1; C2:=0;sol;with(plots): plot(rhs(sol),x=0..10);

                                     C1 := 1

                                     C2 := 0

                                         2           2
                 y(x) = exp(- x) (- 3/2 x  + 1) + 3 x  exp(- x)
# EXAMPLE  #2.
# 
# 
> diffeq:=diff(y(x),x$2)+4*y(x)=tan(x);

                              /   2      \
                              |  d       |
                    diffeq := |----- y(x)| + 4 y(x) = tan(x)
                              |   2      |
                              \ dx       /
> char_eq:=r^2+4; ev:=solve(",r);

                                            2
                                char_eq := r  + 4

                                ev := 2 I, - 2 I
> FS:=[cos(2*x),sin(2*x)]; WW:=Wronskian(FS,x);

                           FS := [cos(2 x), sin(2 x)]

                             [   cos(2 x)     sin(2 x)  ]
                       WW := [                          ]
                             [ - 2 sin(2 x)  2 cos(2 x) ]
> R:=matrix(2,1,[0,tan(x)]);

                                      [    0   ]
                                 R := [        ]
                                      [ tan(x) ]
> (1/WW) &* R; evalm("); map(simplify,");

                                      1
                                    ---- &* R
                                     WW

                         [          sin(2 x) tan(x)    ]
                         [ - 1/2 --------------------- ]
                         [               2           2 ]
                         [       cos(2 x)  + sin(2 x)  ]
                         [                             ]
                         [         cos(2 x) tan(x)     ]
                         [  1/2 ---------------------  ]
                         [              2           2  ]
                         [      cos(2 x)  + sin(2 x)   ]

                         [               2            ]
                         [         cos(x)  - 1        ]
                         [                            ]
                         [              2             ]
                         [     (2 cos(x)  - 1) sin(x) ]
                         [ 1/2 ---------------------- ]
                         [             cos(x)         ]
> map(int,",x)+matrix(2,1,[C1,C2]); U:=evalm(");

                   [    1/2 sin(x) cos(x) - 1/2 x   ]
                   [                                ]   [ 1 ]
                   [             2                  ] + [   ]
                   [ - 1/2 cos(x)  + 1/2 ln(cos(x)) ]   [ 0 ]

                          [  1/2 sin(x) cos(x) - 1/2 x + 1 ]
                          [                                ]
                     U := [             2                  ]
                          [ - 1/2 cos(x)  + 1/2 ln(cos(x)) ]
> sol:=y(x)=dotprod(FS,convert(U,vector)); 

 sol := y(x) = cos(2 x) (1/2 sin(x) cos(x) - 1/2 x + 1) + sin(2 x)

                  2
     (- 1/2 cos(x)  + 1/2 ln(abs(cos(x))) - 1/2 I (1/2 - 1/2 signum(cos(x))) Pi)
> dsolve(diffeq,y(x));

                                                                            2
    y(x) = 1/2 cos(2 x) sin(x) cos(x) - 1/2 cos(2 x) x - 1/2 sin(2 x) cos(x)

         + 1/2 sin(2 x) ln(cos(x)) + _C1 cos(2 x) + _C2 sin(2 x)
# HERE, I used Dsolve to show you how Maple would solve
# the problem, and so you can compare it (and check)
# to see our answer.
# 
> _C1:=1; _C2:=-1;sol; plot(rhs(sol),x=0..10);

                                    _C1 := 1

                                    _C2 := -1

 y(x) = cos(2 x) (1/2 sin(x) cos(x) - 1/2 x + 1) + sin(2 x)

                  2
     (- 1/2 cos(x)  + 1/2 ln(abs(cos(x))) - 1/2 I (1/2 - 1/2 signum(cos(x))) Pi)
>