<PLAINTEXT>
# This demo shows how to use maple to solve higher order linear differential equations.
> restart:
# Example 1. Consider the second order initial value problem given by y" + 2y' + y =x^2 + 1 + e^x, y(0)=0, y'(0)=2.
# The next two steps show how to define the problem by maple.

> diffeq1:=diff(y(x),x$2) + 2*diff(y(x),x) + y(x) = x^2 +1 + exp(x);
> inits1:=y(0)=0, D(y)(0)=2;
# To get the general solution to the differential equation, we use:
> sol1g:=dsolve(diffeq1,y(x));
# Note that the solution above involves two constants, namely _C1 and _C2. We can solve _C1 and _C2 by our initial conditions. 
# To solve the initial value problem, we use:
> dsolve({diffeq1,inits1}, y(x));
# If you don't like the output format above, you may try:

> expand("); sol1p:=simplify(",exp);
# To plot the solution of the initial valued problem on the interval [-1,3], we may use

> plot(rhs(sol1p), x=-1..3);
# To obtain a floating decimal value for the solution at x=2.5, we may use the command "subs" and "evalf":
> subs(x=2.5, sol1p);
> evalf(");
# Example 2. Consider a third order initial valued problem, y''' + y" + 3y'-5y=2+6x-5x^2, y(0)=-1, y'(0)=1, y'''(0)=-3.
# Once again, we first define the differential equation and the initial conditions.
> diffeq2:=diff(y(x),x$3) + diff(y(x), x$2) + 3*diff(y(x),x) - 5*y(x) = 2+6*x-5*x^2;
> inits2:=y(0)=-1, D(y)(0)=1, (D@@2)(y)(0)=-3;
# Then we solve the initial value problem:
> dsolve({diffeq2,inits2},y(x));
# If you prefer a different form of the answer:
> expand(");sol2:=simplify(",exp);
# Example 3. Not every initial value problem can be solved by dsolve. Here is an example that we cannot get the solution by dsolve. 
# We have to approximate the solution numerically.

# First, define the initial value problem:
> diffeq3:=t^2*diff(y(t),t$2)-y(t)*diff(y(t),t)=y(t)^(-3);
> inits3:=y(1)=0.5, D(y)(1)=2.1;
# Solve the initial valued problem numerically, and pay attention to the key word "numeric" in the dsolve command.:
> sol3:=dsolve({diffeq3,inits3},y(t),numeric);
# If you want the approximation value of the solution at x=1.5:
> sol3(1.5);
# If you want to plot the graph of the numerical solution on the interval [0.6, 1.2], then do the following steps:
> with(plots):
> odeplot(sol3,[t,y(t)],0.6..1.2);