Ordinary Differential Equations

 
# FIRST ORDER ORDINARY SCALAR EQUATIONS
# Maple solves differential equations using dsolve:

 de := diff( y(x), x ) = y(x) - sin(x);
 sol:= dsolve( de, y(x));

#  Initial conditions give a unique solution that can be graphed.

 de := diff( y(x), x ) = y(x) - sin(x);
       # Also: f:=(x,y)->y-sin(x); de:=D(y)(x) = f(x,y(x));
       # But: diff(y,x)=y-sin(x) is a SYNTAX ERROR
 ic := y(0)=1;
 sol:= dsolve( [de,ic], y(x)); 
       # Both syntax [de,ic] and syntax {de,ic} will work.

# To graph the solution, extract the expression after the equal sign
# of the answer with Maple function rhs().

 u:=rhs(sol);
 plot(u,x=0..10);

#  Try ?dsolve for more information and examples.

# DIRECTION FIELDS and PHASE PLOTS
# The DEplot command from the DEtools package can plot
# a direction field and threaded curves.
# EXAMPLE. Field only for y'=y-sin(x)

 with(DEtools):
 de := diff(y(x),x) =y(x -sin(x);
 DEplot(de,y(x),x=-3..3,y=-3..3);

# EXAMPLE. Direction field and threaded curves for y'=y-sin(x)

 ics:=[ [y(0)=0], [y(0)=1], [y(0)=3] ];
   # Set of initial conditions. A double list.
   # Identical Result: ics:=[[0,0],[0,1],[0,3]]:
 opts:=color='black',linecolor='red',arrows='line',dirgrid=[30,30]:
 DEplot(de,y(x),x=-3..3,y=-3..3,opts,ics);

# DSOLVE INTERACTIVE TOOL
# To find the tool, start maple, then goto menu item
#   TOOLS ==>  ASSISTANTS ==> ODE ANALYZER
# Once the tool starts, enter the differential equation as
##   diff(x(t),t)=x(t)*(1-x(t))
# Fill in initial condition x(0)=1.2 as two numbers
##   0 1.2
# Try buttons 
   ## SOLVE SYMBOLICALLY
     ##  Click on SHOW MAPLE COMMANDS
     ##  Click on SOLVE
     ##  Click on PLOT 

# SECOND ORDER ORDINARY SCALAR EQUATIONS

# Solve y'' + y = 0 with initial conditions

  de1 := diff( y(x),x,x) + y(x) = 0;
  ic := y(0)=1,D(y)(0)=0;
  dsolve( [de1,ic], y(x) );

# Solve y'' + y = sin(x) with initial conditions

  de2 := diff( y(x),x,x) + y(x) = sin(x);
  ic := y(0)=1,D(y)(0)=0;
  dsolve( [de2,ic], y(x) );

# THIRD ORDER ORDINARY SCALAR EQUATIONS

# Solve y''' + 2y'' + y' = x^2 without initial conditions
# Syntax   D@D(y)(x) is the same as   diff(y(x),x,x)
#        D@D@D(y)(x) is the same as diff(y(x),x,x,x)

  L:=D@D@D+2*D@D+D; de3:=L(y)(x)=x^2;
  dsolve(de3,y(x));

# DYNAMICAL SYSTEMS 

# Solve x'=2x+3y, y'=4y, x(0)=1, y(0)=2
# Method: Use matrices and DEtools
 A:=Matrix([[2,3],[0,4]]); # Matrix of coefficients
 sol:=DEtools[matrixDE](A,t); # Solve u'(t)=Au(t)
 Phi:=unapply(sol[1],t);      # Fundamental matrix Phi(t)
 u:=Phi(t).Phi(0)^(-1).<1,2>; # Solution of u'=Au, u(0)=<1,2>

# Solve x'=2x+3y + exp(t), y'=4y, x(0)=1, y(0)=2
# Method: dsolve. See also DEtools[matrixDE].
 A:=Matrix([[2,3],[0,4]]);   # Matrix of coefficients
 F:=< exp(t),0>;             # Define F in u'=Au+F
 u:=unapply(< x(t),y(t)>,t); # Symbolic solution variables
 sys:=map(diff,u(t),t)-A.u(t)-F;  # DEs are sys[1]=0, sys[2]=0
 ic:=u(0)-<1,2>;             # initial conditions ic[1]=0, ic[2]=0
 dsolve({sys[1],sys[2],ic[1],ic[2]},{x(t),y(t)});
 # Ans: x(t) = 3*exp(4*t)-exp(t)-exp(2*t), y(t) = 2*exp(4*t)

# Find exp(At)

 A:=Matrix([[2,3],[0,4]]); 
 LinearAlgebra[MatrixExponential](A,t);
  

Partial Differential Equations

EXAMPLE: Heat Equation, Numeric

 
# Solve a partial differential equation numerically.

# Heat Equation u_t = (0.1) u_{xx}

 pde:=diff(u(x, t), t)=(1/10)*(diff(u(x, t), x, x));

# Initial boundary conditions
# Math boundary conditions: u=1 at t=0,  u=0 at x=0,  u_x=0 at x=1
# Symbols: D[1] is partial on x, D[2] is partial on t

 ics:= u(x, 0) = 1, u(0, t) = 0, (D[1](u))(1, t) = 0;

# Solve the equation numerically,

 sol := pdsolve(pde,{ics},numeric);

# Plot temperature snapshots on one graphic

 p1 := sol[plot](t = 1/10, color = black):
 p2 := sol[plot](t = 1/5, color = red):
 p3 := sol[plot](t = 3/4, color = blue):
 opts:=title = "Temperature snapshots for t = 0.1, 0.5, 0.75":
 plots[display]([p1, p2, p3], opts);

Reducing PDE Display Complexity

# with(PDEtools);
 PDE:=diff(u(x,t)*v(t)*w(y),x)+diff(u(x,t)*v(t)*w(y),y,t);
 Vars := [u(x,t),v(t),w(y)];
 PDEtools[ToJet](PDE,Vars,jetnotation=jetvariables);
# The PDE will be re-written in the commonly used notation
#           u[x]*v*w+u[t]*v*w[y]+u*v[t]*w[y]
# with pretty-printing using subscripts, e.g., u[x]-->u_x