# CHAPTER 8 - DIFFERENTIAL EQUATIONS
#
# 8.1   Solving ordinary differential equations
#
#
#
# MAPLE SESSION 1
#
>   y:='y': 
>   Y:=y(x);
>   dY := diff(Y,x);
>   ddY := diff(%,x);
>   de := ddY+5*dY+6*Y = sin(x)*exp(-3*x); 
>   ans := dsolve(de, Y);
#
#
# MAPLE SESSION 2
#
>    	sol := rhs(ans);
>    	dsol:=diff(sol,x);
>    	ddsol:=diff(sol,x,x);
>    	ddsol + 5*dsol + 6*sol;
#
# 8.1.1   Implicit solutions
#
#
#
# MAPLE SESSION 3
#
>    DE:=(3*y(x)^2+exp(x))*diff(y(x),x)+exp(x)*(y(x)+1)
>  +cos(x) = 0;
>    dsolve(DE,y(x),implicit);
#
#
# MAPLE SESSION 4
#
>    dsolve(DE,y(x));
#
# 8.1.2   Initial conditions
#
#
#
# MAPLE SESSION 5
#
>       de;
>       dsolve({de,y(0)=-5/2,D(y)(0)=2},Y);
#
#
# MAPLE SESSION 6
#
>    DE := y(x) +D(y)(x) -3*(D\@\@2)(y)(x)+ (D\@\@3)(y)(x) 
>  = exp(-x)*(10-4*x);
>    dsolve({DE,y(0)=5,D(y)(0)=6,(D\@\@2)(y)(0)=3},y(x));
#
# 8.1.3   Systems of differential equations
#
#
#
# MAPLE SESSION 7
#
>    de1 := diff(y(x),x) + diff(z(x),x) = exp(x); 
>    de2 := diff(y(x),x) -3*z(x) = x; 
>    dsolve({de1,de2,y(0)=8/9,z(0)=10/9},{y(x),z(x)}); 
#
# 8.2   First order differential equations
# 8.2.1   {\tt odeadvisor}
#
#
#
# MAPLE SESSION 8
#
>    restart: 
>    with(DEtools): 
>    DE := diff(y(x),x) = (x-3)/y(x)^2;
>    odeadvisor(DE);
#
#
# MAPLE SESSION 9
#
>   odeadvisor(DE,help);
#
#
# MAPLE SESSION 10
#
>    dsolve(DE, y(x)); 
#
#
# MAPLE SESSION 11
#
>    with(DEtools): 
>    DE := diff(y(x),x) = (x-3)/exp(y(x));
>    odeadvisor(DE);
>    infolevel[dsolve]:=3;
>    dsolve(DE,y(x));
#
# 8.2.2   Integrating factors
#
#
#
# MAPLE SESSION 12
#
>    with(DEtools): 
>    de:=(x^4+y(x))+x*(y(x)^2+ln(x))*diff(y(x),x)=0;
>    mu := intfactor(de,y(x));
#
#
# MAPLE SESSION 13
#
>    de2:=expand(mu*de);
>    odeadvisor(de2);
#
# 8.2.3   Direction fields
#
#
#
# MAPLE SESSION 14
#
>    restart: 
>    with(DEtools): 
>    DE:= diff(y(x),x)=4*x/y(x);
>    DEplot(DE,y(x),x=-2..2,y=-1..3); 
#
#
# MAPLE SESSION 15
#
>    DEplot(DE,y(x),x=-2..2,[[y(0)=2]],y=-1..3); 
#
#
# MAPLE SESSION 16
#
>    DEplot(DE,y(x),x=-2..2,[[y(0)=2],[y(1)=1]],y=-1..3); 
#
#
# MAPLE SESSION 17
#
>    DEplot(DE,y(x),x=-2..2,[[y(0)=2],[y(1)=1]],y=-1..3,
>  stepsize=0.001); 
#
#
# MAPLE SESSION 18
#
>    DEplot(DE,y(x),x=-2..2,[[y(0)=2],[y(1)=1]],y=-1..3,
>  stepsize=0.001, arrows=NONE);
#
# 8.3   Numerical solutions
#
#
#
# MAPLE SESSION 19
#
>    de:=diff(y(x),x)=x^2+y(x)^3: 
>    IVP:={de,y(0)=0};
>    nsol := dsolve(IVP,y(x),numeric);
#
#
# MAPLE SESSION 20
#
>    for i from 0 to 9 do 
>      nsol(i/10.);
>    end do;
#
#
# MAPLE SESSION 21
#
>   with(plots): 
>   odeplot(nsol,[x,y(x)],0..1,title="Numerical solution to 
>   dy/dx=x^2+y^3, y(0)=0"); 
#
# 8.4   Second and higher order linear de's
# 8.4.1   Constant coefficients
#
#
#
# MAPLE SESSION 22
#
>    with(DEtools): 
>    de:=diff(y(x),x,x)+2*diff(y(x),x)+4*y(x)=0;
>    constcoeffsols(de,y(x));
#
# 8.4.2   Variation of parameters
#
#
#
# MAPLE SESSION 23
#
>    with(DEtools): 
>    varparam([cos(x),sin(x)],tan(x),x);
#
# 8.4.3   Reduction of order
#
#
#
# MAPLE SESSION 24
#
>     with(DEtools): 
>     de:=diff(y(x),x,x,x)+diff(y(x),x,x)+3*diff(y(x),x)
>  -5*y(x)=0;
>     reduceOrder(de,y(x),exp(x));                                
>     dsolve(%,y(x));
#
# 8.5   Series solutions
#
#
#
# MAPLE SESSION 25
#
>    de:=(1+x^2)*diff(y(x),x,x) + 3*x*diff(y(x),x) + y(x) 
>  = 0;
>    dsolve(de,y(x),type=series);
#
#
# MAPLE SESSION 26
#
>    dsolve({de,y(0)=1,D(y)(0)=0},y(x),type=series);
#
#
# MAPLE SESSION 27
#
>    Order:=10;
>    dsolve({de,y(0)=1,D(y)(0)=0},y(x),type=series);
#
#
# MAPLE SESSION 28
#
>   Order := 6:
>    de:=(x^2-3*x)*diff(y(x),x,x) + 2*x*diff(y(x),x) + y(x) 
>  = 0;
>    dsol := dsolve(de,y(1)=4,D(y)(1)=-2,y(x),type=series);
#
#
# MAPLE SESSION 29
#
>    psol:=convert(rhs(dsol), polynom);
>    plot(psol,x=0..2);
#
# 8.5.1   The method of Frobenius
#
#
#
# MAPLE SESSION 30
#
>     with(DEtools): 
>     DE:=t*(1-t)*diff(y(t),t,t)+(c - (1+a+b)*t)*diff(y(t),t)
>  -a*b*y(t)=0;
>   regularsp(DE,t,y(t));
>     indicialeq(DE,t,0,y(t));
>     solve(%,t);
#
# 8.6   The Laplace transform
#
#
#
# MAPLE SESSION 31
#
>    with(inttrans);
#
#
# MAPLE SESSION 32
#
>    laplace(t*exp(3*t),t,s);
#
#
# MAPLE SESSION 33
#
>    invlaplace((s-1)/(s^2+s+6),s,t);
#
# 8.6.1   The Heaviside function
#
#
#
# MAPLE SESSION 34
#
>    plot(Heaviside(t),t=-1..1,discont=true,thickness=3,
>  title="The Heaviside Function H(t)"); 
#
#
# MAPLE SESSION 35
#
>     restart: 
>     with(inttrans): 
>     dy:=diff(y(t),t);
>     ddy:=diff(dy,t);
>     dddy:=diff(ddy,t);
>     addtable(laplace,y(t),Y(s),t,s);
>     laplace(y(t),t,s);
>     laplace(dy,t,s);
>     laplace(ddy,t,s);
#
#
# MAPLE SESSION 36
#
>   f := sin(t)*(Heaviside(t-Pi)-Heaviside(t-2*Pi));
#
#
# MAPLE SESSION 37
#
>   y(0) := 0;
>   D(y)(0) := 0;
#
#
# MAPLE SESSION 38
#
>   ddy + 4*y(t);
#
#
# MAPLE SESSION 39
#
>   laplace(ddy+4*y(t),t,s) = laplace(f,t,s);
>     simplify(%);
#
#
# MAPLE SESSION 40
#
>     solve(%,Y(s));
>     S := factor(%);
#
#
# MAPLE SESSION 41
#
>     ysol := invlaplace(S,s,t);
#
#
# MAPLE SESSION 42
#
>     y:='y': D(y):='D(y)':
>     IC := y(0) = 0, D(y)(0) = 0;
>     DE := ddy + 4*y(t) = f:
>     dsolve({DE,IC},y(t));
>     ysol2:=rhs(%);
>     simplify(ysol - ysol2);
#
# 8.6.2   The Dirac delta function
#
#
#
# MAPLE SESSION 43
#
>    with(inttrans): 
>    assume(t0>=0): 
>    laplace(Dirac(t-t0),t,s);
#
#
# MAPLE SESSION 44
#
>     restart: 
>     with(inttrans): 
>     dy:=diff(y(t),t): 
>     ddy:=diff(dy,t): 
>     addtable(laplace,y(t),Y(s),t,s): 
>     f := exp(-t) + 2*Dirac(t-1): 
>     y(0) := 0: 
>     D(y)(0) := -1: 
>     LS := ddy + 2*dy + y(t);
#
#
# MAPLE SESSION 45
#
>     laplace(LS,t,s) = laplace(f,t,s);
#
#
# MAPLE SESSION 46
#
>     solve(%,Y(s));
>     S := factor(%);
#
#
# MAPLE SESSION 47
#
>     invlaplace(S,s,t);