# Several DEs from Boyce and Diprima are solved and checked 
# using Maple.  
# The use of DEplot1 is introduced as well.
#
# Art Belmonte
# Mon, 27/May/96
# Math 308-509 [Maple V Release 4] 
# Examples
# 
# Ex. 1 (S2.1, T-22/1): Recall that `T' stands for your textbook (`M'
# would signify lab manual), 22/1 means page 22 / problem #1, and that
# S2.1 indicates the problem is from Section 2.1. Note the arbitrary
# multiplicative constant _C1 (we used just C when doing the problem by
# hand). This indicates the differential equation has an infinite number
# of solutions (see 22/1 revisited below).
# 
# The naming of DE(s) and/or IC(s) permits future reference. The
# "unapplication" checking technique may appear cryptic at first, but
# allows for uniform and simultaneous handling of checks on DEs and ICs
# (see other examples below). Both of these techniques are from the new
# Maple V Release 4 _Learning Guide_. Wrap the deq (or IC or IVP, in
# other problems) in the check with simplify if necessary. Your mileage
# may vary...
#
> with(DEtools):
#
> unassign('y');
> deq:=diff(y(x), x) + 3*y(x) = x + exp(-2*x);

                     /d      \
              deq := |-- y(x)| + 3 y(x) = x + exp(-2 x)
                     \dx     /

> sol:=dsolve(deq, y(x));

        sol := y(x) = 1/3 x - 1/9 + exp(-2 x) + exp(-3 x) _C1

> y:=unapply(subs(sol, y(x)), x);

          y := x -> 1/3 x - 1/9 + exp(-2 x) + exp(-3 x) _C1

> check:=deq;

                check := x + exp(-2 x) = x + exp(-2 x)

# T-22/1 revisited! This illustrates the direction field for a
# differential equation (with or without integral solution curves). In
# order to use the DEplot command, the DEtools package must be loaded.
# This was done via the read command at the beginning of the worksheet.
> unassign('y');
> DEplot(deq, y(x), x=-2..2, y=-2..2);
> inits:=[ [0, 0], [0, 1], [0, -1] ];

                  inits := [[0, 0], [0, 1], [0, -1]]

> DEplot(deq, y(x), x=-2..2, y=-2..2, inits);
# Ex. 2 (S2.2, T-28/5)
> unassign('y');
> deq:=x*diff(y(x), x) + 2*y(x) = x^2 - x + 1;

                        /d      \             2
               deq := x |-- y(x)| + 2 y(x) = x  - x + 1
                        \dx     /

> IC:=y(1) = 1/2;

                           IC := y(1) = 1/2

> IVP:={deq, IC};

                              /d      \             2
        IVP := {y(1) = 1/2, x |-- y(x)| + 2 y(x) = x  - x + 1}
                              \dx     /

> sol:=dsolve(IVP, y(x));

                                    4      3      2
                                 3 x  - 4 x  + 6 x  + 1
              sol := y(x) = 1/12 ----------------------
                                            2
                                           x

> y:=unapply(subs(sol, y(x)), x);

                                  4      3      2
                               3 x  - 4 x  + 6 x  + 1
                y := x -> 1/12 ----------------------
                                          2
                                         x

> check:=simplify(IVP);

                                  2            2
            check := {1/2 = 1/2, x  - x + 1 = x  - x + 1}

# Ex. 3 (S2.3, T-36/9): Note that only ONE of the solutions solves the
# IVP, y(0)=1. Be careful and check your work!
> unassign('y');
> deq:=x + y(x)*exp(-x)*diff(y(x), x) = 0;

                                        /d      \
                deq := x + y(x) exp(-x) |-- y(x)| = 0
                                        \dx     /

> IC:=y(0) = 1;

                            IC := y(0) = 1

> IVP:={deq, IC};

                                   /d      \
          IVP := {x + y(x) exp(-x) |-- y(x)| = 0, y(0) = 1}
                                   \dx     /

> sol:=dsolve(IVP, y(x));

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

> y:=unapply(subs(sol, y(x)), x);

                                                    1/2
              y := x -> (-2 x exp(x) + 2 exp(x) - 1)

> check:=simplify(IVP);

                       check := {1 = 1, 0 = 0}

# T: 36/9 --- Approximate interval in which the solution is valid.
> expr:=rhs(sol)^2; need:=expr > 0;

                  expr := -2 x exp(x) + 2 exp(x) - 1


                need := 0 < -2 x exp(x) + 2 exp(x) - 1

> plot(expr, x=-3..1);
> left_endpt:=fsolve(expr=0, x=-2..-1);

                      left_endpt := -1.678346990

> right_endpt:=fsolve(expr=0, x=0..1);

                      right_endpt := .7680390470

>