Chapter 1
Introduction
1.1 Motivational example s
To solve real-world problems, scientists formulate them in terms of mathematical expressions.
Those expressions are usually in the form of either ordinary or partial differential equa-
tions (DEs). This book discusses the first type of eq uations, ordinary differential equations
(ODEs).
1.1.1 Population g rowth
In 1798, T. Malthus presented a mathematical mod el for the growth of population.
According to his model, the population increases exponentially that follows the following
equation
dP
dt
= rP :
In other words, the rate of change of the population function P (t) at time t th a t is
dP
dt
(t)
is proportional to the present population P (t), and the constant of the proportionality is a
constant r > 0 that expresses the off-spring or reproduction rate. It is simply seen that P (t)
is exponential of the form P (t) = P (0) e
rt
. On the other hand, fo od production (acco rding
to agricultural development) increases linearly, and thus, one leads to a pessimistic view of
the future of humankind in starvation and misery. This was a shock at that time. In 1831,
J. P. Verhulst published a paper and showed that the Malthus mathematical model is
not accurate. He proposed a new model for
dP
dt
as follows
dP
dt
= rP
1 ¡
P
K
:
(1.1)
where K > 0 is called carrying capacity and reflects the limitations of the population growth
imposed by the environment. It is simply verified by the direct substitution that P (t) is the
following function if it follows the above differential equation
P (t) =
KP
0
P
0
+ (K ¡P
0
)e
¡rt
;
1
where P
0
= P (0) is the initial population at time t = 0.
1.1.2 A falling body
Consider a mass m falling from a height H. Assume that f
r
, the air resistance against the
bo dy is proportional to its velocity f
r
= ¡kv, where k > 0 is a constant (the negative sign
enters because the resistance always acts against the direction of motion)
f
r
= ¡kv(t)
m
f = mg
We would like to determine h, the height of the mass at any instance of time t. Applying
the Newton's second law m
dv
dt
= f results to the following equation for v(t), the velocity
of the mass
m
dv
dt
= mg ¡kv: (1.2)
The above equation is an ordinary differential equation that describes the rate of change of
v(t) in terms of t as well as v(t) itself. If we solve the above equation and derive v(t), then
the h(t), the height of the mass at a ny instance of time is simply determined by the following
integral
h(t) = h
0
¡
Z
0
t
v(s) ds;
where h
0
= h(0) is th e initial height of the m a ss. It is simply verified that the following
function solves Eq.1.2 for arbitrary parameter c
v(t) = ce
¡
kt
m
+
mg
k
:
In fact, by taking derivative, we have
dv
dt
=¡
ck
m
e
¡
kt
m
, and substituting this into the differential
equation, transforms the equation to an identity in terms of t. The obtained solutio n that
contains a parameter c i s called the general solution of the equation. The parameter c is
uniquely determined if an initial condition v(0) is prov ided for the equation. For example,
if v(0) = 0, then c = ¡
mg
k
. For this value of c, the solution becomes
v(t) =
mg
k
1 ¡e
¡
kt
m
:
2 Introduction
This is called a particular solution to the equation. If the initial condition v(0) is set to v
0
,
then
v(0) = c +
mg
k
= v
0
;
and then the solution becomes
v(t) =
v
0
¡
mg
k
e
¡
kt
m
+
mg
k
:
Problem 1.1. Solve th e above differential equation by integration as follows
Z
v
0
v
ds
¡
k
m
s + g
=
Z
0
t
ds;
where v
0
= v(0) the initial velocity of the mass.
Problem 1.2. For high speed motions, the drag force acts as f
r
= ¡kv
2
instead of linear term ¡kv.
Assume a particle under the constant force f
0
is experiencing the drag force of form ¡kv
2
. The differential
equation in this case is a s follows
m
dv
dt
= ¡kv
2
+ f
0
:
Rewrite the equation as
dv
¡
k
m
v
2
+
f
0
m
= dt;
and integrate bo th sides of the equation to find the general solution v(t). If v(0) = 0, find the particular
solution.
1.1.3 The geo metry of the equation
Consider again the equation (1.2) and rewrite it as
dv
dt
= g ¡
k
m
v:
The physics of th e equation states that the rate of change of the velocity of the falling body
at a ny instance of time,
dv
dt
(t), is equal to g ¡
k
m
v(t). For example, at time t = 0, the change
of the velocity is equal to
dv
dt
(0) = g ¡
k
m
v(0);
and if v(0) = v
0
given , then
dv
dt
(0) = g ¡
k
m
v
0
.
Remembe r that
dv
dt
(t) represents geometrically the slope of tangent vector to the function
v =v(t) at time t. The following figure shows two solution curves of the equation for v(0)= 10
and v(0) = 25 for m = 1; k = 0.5. The tangent vectors are shown by arrows in the plane (t; v).
The curves are tangent everywhere to the vectors according to the differential equation
v
0
= g ¡0.5v.
1.1 Motivational examples 3
0 2 4 6 8 10
10
12
14
16
18
20
22
24
26
28
30
v
0
=25
v
0
=10
Observe that both solutions tends to the same value 2g, that is,
lim
t!1
v(t) = 20:
Problem 1.3. Use a computer software and draw the slope field of equation
dv
dt
= ¡kv
2
+ f
0
for your
choices of k > 0, a nd f
0
> 0. Find the steady state so lution to the equation.
1.1.4 A variable mass
Consider a raindrop in the shape of a ball falling straight down, and suppose it absorbs water
in th e form of the stem from the ambient space. Let us obtain the equation of the motion of
such a raindrop. Assume the raindrop is in the shape of a ball with the initial mass m
0
and
that its absorption rate is proportional to its surface area. In this ca se, the mass changes
with time and the Newton's law reads as follows
d
dt
(mv) = ¡mg; (1.3)
where m = m(t) changes with time. According to th e assumption on the shape of the rain
drop, and the fo rmula S = 4πr
2
, one can w rite
dm
dt
= kr
2
;
for some k > 0 . Note that the mass of a ball of vo lume V and density ρ is m = ρV =
4πρ
3
r
3
,
and therefore r
2
=
3
4πρ
m
2
3
, that by the assumption
dm
dt
= kr
2
gives
dm
dt
= k
3
4πρ
2
3
m
2
3
:
4 Introduction
Therefore, we obtain the following system of equations
8
>
>
<
>
>
:
m
dv
dt
+ v
dm
dt
= ¡mg
dm
dt
= k
3
4πρ
2
3
m
2
3
:
However, the equation (1.3) can be integrated as
Z
0
t
d(mv) = ¡g
Z
0
t
m(s) ds;
that yields
m(t) v(t) ¡m(0)v(0) = ¡g
Z
0
t
m(s) ds:
Here we n eed m(t) to calculate the right hand side of the above equality. We can rewrite
the equation for m(t) as follows
m
¡
2
3
dm = k
3
4πρ
2
3
dt;
and therefore
d
¡
3m
1
3
= k
3
4πρ
2
3
dt;
that is solved by integration as
3
h
m(t)
1
3
¡m(0)
1
3
i
= k
3
4πρ
2
3
t;
or
m(t) =
k
1
t + m
0
1
3
3
;
where k
1
=
k
3
3
4πρ
2
3
. If v(0), the initial velocity of the rain drop is zero, we o btain
mv = ¡g
Z
0
t
k
1
s + m
0
1
3
3
ds = ¡
g
4k
1
h
k
1
t + m
0
1
3
4
¡m
0
4
3
i
;
and finally
v(t) = ¡
g
4k
1
"
m ¡
m
0
4
3
m
#
:
Problem 1.4. The second Newton law reads m
d
2
x
dt
2
= f only if m is constant. When m is variable, the
formula needs some modification. First, we write the Newton's second law in its original form
dp
dt
= f ,
where p is the momentum of the system.
i. Consider a system consisting two m asses: the mass m that moves with velocity v and mass δm
that moves with velocity u. The momentum of the system is
p(t) = mv(t) + δmu(t): (1.4)
Suppose that at time t + δt, the mass δm att aches to mass m and the combined mass moves with
velocity v(t) + δv. Find
dp
dt
(t) by the following limit
dp
dt
(t) = lim
δt! 0
p(t + δt) ¡ p (t)
δt
:
ii. Conclude that the Newton's second law in this case read as follows
d
dt
(mv) ¡u
dm
dt
= f :
1.1 Motivational examples 5
iii. Use the above formula and derive
d
dt
(mv) = ¡mg we used in the above example.
1.1.5 A geometrical problem
Example 1.1. We are going to make a mirror that collects light rays coming from a far
source into a focal point.
focal
φ
x
φ
φ
x
y
2φ
y
To obtain the shap e of the mirror, we have to make a differential equation as follows. Let
y(x) be the shape function of the mirror. By the above figure, we can w rite
dy
dx
= tanφ , a nd
since tan(2φ) =
y
x
, we obtain by the trig o nometric identity
tan(2φ) =
2 tanφ
1 ¡tan
2
φ
=
2y
0
1 ¡y
0
2
;
and finally we re a ch the following equation fo r the shape of the mirror
(1 ¡y
0
2
)y = 2xy
0
:
The above equation is called the Lagrange equation. It is simply seen that function
x =
y
2
4c
¡c;
solves the equation for any c =/ 0. Let u s verify the claim. Taking derivative of the function
with respect to x g ives 1 =
1
2c
yy
0
, and thus y
0
=
2c
y
. Substituting this into the differential
equation results to y
2
¡4c
2
= 4 cx , and substituting x from the solution we get an identity.
Therefore, the obtained solution transforms the differential equation into an identity. Note
that the shape of the mirror is a parabola.
Problem 1.5. Let us try to solve the above differential equation. If we denote y
0
by p, then y is
represented by the parametric function y =
2xp
1 ¡ p
2
. We need to find x = x(p).
i. Find a differential equation for x in terms of p by calculating
dx
dp
. You should find it as follows
p
dx
dp
=
¡2
1 ¡ p
2
x:
ii. Solve the differential equa tion obtain in part (i) by integrating
Z
dx
x
=
Z
¡2
p(1 ¡ p
2
)
dp
6 Introduction
and find x(p)
iii. Eliminate parameter p from system x(p) and y =
2xp
1 ¡ p
2
and derive the solution.
1.2 Preli minaries
1. (Differential) The differential of a single variable function y = f(x) is defined by
the relation dy = f
0
(x) dx, where dx is the differential of x and f
0
(x) is the derivative
function of f which is defi ned as follows
f
0
(x) = lim
h!0
f(x + h) ¡f(x)
h
:
For example, if y = e
αx
, then dy = αe
αx
dx. The fundamental theorem of calculus states
a formula for the integration of d ifferential, that is,
Z
x
0
x
dy =
Z
x
0
x
f(x) dx = f(x) ¡ f(x
0
):
2. (Tota l differential) For a two-variable smooth functions z = f (x; y), the differential
dz is defined as
dz =
@f
@x
dx +
@f
@y
dy;
where
@f
@x
;
@f
@y
are partial derivative fu nctions of f
@f
@x
(x; y) = lim
h!0
f(x + h; y) ¡f(x; y)
h
;
@f
@y
(x; y) = lim
h!0
f(x; y + h) ¡f(x; y)
h
:
For example, if f (x; y) = x
2
+ y
2
+ sin(xy) then
df = (2x + y cos(xy)) dx + (2y + x cos(xy))dy:
The fundamental theorem of calculus or the integration over a path that connects
points (x
0
; y
0
) and (x; y) is
Z
(x
0
;y
0
)
(x;y)
dz = f (x; y ) ¡ f (x
0
; y
0
):
3. (Classical differentia l equation) In the period of the early development of Cal-
culus by I. Newton in seventeen century, mathematicians studied equations of the
following form
M(x; y) dx + N(x; y) dy = 0; (1.5)
1.2 Preliminaries 7
for given functions M ; N . These types of equations called differential equations show
themselves in the process of the mathematical modeling of real-world problems, in
physics, biology, chemistry, engineering, and geometry. The goal is to find a contin-
uously differential function y = y(x) in an open interval I such that
M(x; y(x)) dx + N(x; y(x)) y
0
(x) dx = 0;
or equivale ntly
M(x; y(x)) + N(x; y(x)) y
0
(x) = 0
for all x 2I. In several cases, the function y(x) is obtained in the implicit form as
f(x; y) = const:
For example, the equation
df = (2x + y cos(xy)) dx + (2y + x cos(xy))dy = 0;
has the implicit solution
f(x; y) = x
2
+ y
2
+ sin(xy) = c:
5
5
5
5
10
10
10
10
10
15
15
15
15
15
15
15
20
20
20
20
25
25
25
25
-4 -3 -2 -1 0 1 2 3 4
-4
-3
-2
-1
0
1
2
3
4
4. (Implicit functi on theorem) Remember that implicit functions may represents a
true explicit functions y = y(x) or not depending on the implicit function theorem.
Theorem 1.1. Given an implicit function f (x; y) = 0, assume
i. there is a point (x
0
; y
0
) such that f(x
0
; y
0
) = 0,
ii. the partial derivative
@f
@y
is a continuous at (x
0
; y
0
) and
@f
@y
(x
0
; y
0
) =/ 0,
then there is an open interval I = (x
0
¡δ; x
0
+ δ) and a function y: I !R such that
f(x; y(x)) = 0 for all x 2I.
8 Introduction
For example, consider the function f := x
2
+ y
2
+ sin(xy) = 1. At (x
0
; y) = (1; 0),
we have
@f
@y
(1; 0) = 2x + x cos(xy)j
(1;0)
= 3;
and thus there is an interval I = (1 ¡δ; 1 + δ) and a continuo us function y: I !R such
that
x
2
+ y(x)
2
+ sin(xy(x)) = 1; 8x 2I:
5. (Relation with physics) Consid er a force field F =
f(x; y)
g(x; y)
and a mass m moving
under the influence of this field. The path or trajectory of the mass in the (x; y)-plane
is usually represented by a parametric curve γ(t), where t denote the time variable.
The total energy of the mass is expressed by the following expression
E =
1
2
m jvj
2
+ V (x; y);
where V is the potential of the force field
f = ¡
@V
@x
; g = ¡
@V
@y
;
and v is the velocity of the mass. It is simply seen that the total energy of the mass
along its trajectory γ(t) is constant. In fact, we have
dE
dt
= mv ·v
0
+
@V
@x
dx
dt
+
@V
@y
dy
dt
= mv ·v
0
¡f
dx
dt
¡g
dy
dt
= mv ·v
0
¡F ·v = v ·(mv
0
¡F );
where the last term is the immediate result of the chain rule
d
dt
V (x(t); y(t)) and the
observation that the velocity vector v is equal to
x
0
(t)
y
0
(t)
!
. According to the second
Newton's law mv
0
= F , we obtain
dE
dt
= 0. On the other hand, we have
dE
dt
=
@E
@x
dx
dt
+
@E
@x
dx
dt
;
and thus we conclude
dE :=
@E
@x
dx +
@E
@x
dy = 0;
in the (x; y)-plane.
6. (General form of first-order equation s) It turns out in a while that there are
several real-world problems whose differential equation can not be expressed in terms
of the above form studied by Newton and others. For example, equation y =
2xy
0
1 ¡ y
0
2
that we s tud ied in the previous section is of this form. For this reason, mathemati-
cians usually defines the general form of a first-order ordinary differential equation as
follows
G(x; y; y
0
) = 0; (1.6)
where y
0
:=
dy
dx
, a nd G is a functional relation between x; y and y
0
. If the functional
relation can be solved algebraically for y
0
as the explicit form y
0
= f ( x; y), then it can
be put in th e differential form dy ¡ f (x; y)dx = 0.
1.2 Preliminaries 9
7. (Initial valu e problem) An initial value problem associated to a general first-order
differential equation is as follows
G(x; y; y
0
) = 0
y(x
0
) = y
0
;
where (x
0
; y
0
) is called the initial condition for the solution of the equation. For
example, consider the following equation that states the relation between the rate of
change of the population P (t) of a living species and the population itself
dP
dt
= rP
1 ¡
P
K
;
where r is the off-spring rate of the species, and K is the carrying capacity of the
environment or the ambient space. It is simply verified that any function of the form
P (t) =
cK
c + e
¡rt
;
solves the equation for arbitrary parameter c. H owever, if we know the initial pop-
ulation of the species at time t = 0, say P (0) = P
0
, then the obtained solution must
satisfies this initial condition as well. If we put P = P
0
for t = 0, we reach
cK
1 + c
= P
0
;
or equivale ntly, c =
P
0
K ¡ P
0
, and thus the particular solution is
P (t) =
P
0
K
K ¡ P
0
P
0
K ¡ P
0
+ e
¡rt
:
Problem 1.6. What happen if P
0
= K? State the physical implication of this condition.
8. (System of equations) While the equation y
0
= f(t; y) states the rel a tion between
the rate of change of a quantity y in terms of the quantity y and time t, the real-world
problems involve several factors that cause the change in a quantity. For example, the
change of the p o pulation of a certain living species, say prey, depends in addition to its
off-spring rate, the population of other species, say predator. A simple mathematical
model that express the dynamics is as follows
8
<
:
dx
dt
= r
1
x ¡αx
2
¡k
1
xy
dy
dt
= ¡r
2
y + k
2
xy
;
where r
1
; r
2
; α; k
1
; k
2
are positive constants, x(t) is the population of prey, say rabbits,
and y(t) is the population of predator, say foxes. Observe that
d
dt
x(t) depends on y(t)
in addition to x(t), and similarly
d
dt
y(t). For this reason, one cal above equations, a
system of differential equations.
10 Introduction
Problems
Problem 1.7. Let y(x) = xe
x
+ 1. Find dy as a function of dx and determine d y if dx = 3.
Problem 1.8. Consider a disk of radius r. If the rate of change of r is as follows
dr
dt
= R ¡r;
for some R > 0, find lim
t! 1
A(t) where A is the area of the disk.
Problem 1.9. Consider the falling body problem discussed in this section and assume that the mass is
unit, m = 1, and the air resistance is in the form f
r
= ¡kv
2
(in the case when speed is too high, the air
resistance applies in the square of speed).
a) Write down the initial value problem for v assuming v(0) = 0.
b) Solve the equation and show that the solution is
v(t) =
g
k
r
e
2 kg
p
t
¡1
e
2 kg
p
t
+ 1
:
c) Find the terminal velocity of the body if the mass is released from a very high position
Problem 1.10. Consider a water tank in the shape of a cube with dimensions (L; W ; H)
l = H
W
L
H
Q
i
Let h(t) denote the water level in the tank at time t. We would like to determine the water level h(t)
for t > 0 if a current of water run into the tank with the rate Q(t) = k (H ¡h(t))m
3
/s, assumin g h(0) = 0.
a) If V (t) denotes the volume of water in the tank at time t, then the conservation principle states
dV
dt
= input rate ¡ou tpu t rate:
By the above relation, derive a the following differential equation for h(t)
dh
dt
(t) =
k (H ¡h( t))
LW
:
b) Solve the equation by integration and conclude that the water tank does not become full in any
finite time.
Problem 1.11. Consider the water tank problem discussed in the previous problem. Assum e that a
hole is placed at the bottom of the tank that let water runs out with the rate Q
o
= α h
p
where α > 0 is
a constant (We will see later on that this assumption makes sense due to the Torricelli law).
a) Write down the differential equation describing the water level h(t) in the tank.
b) The terminal level or equilibrium is defin ed as the water level reaches to a certain level without
any further change. Find the termi na l level of the water in the tank and determine whether the
tank become full or not.
Problem 1.12. In the water tank problem 1.10, assume L = W = 2, H = 1 and Q
i
= sin
2
φ.
i. Write the equation for the water level h(t) and solve it to find h(t) for t > 0 assuming h(0) = 0.
1.2 Preliminaries 11
ii. Find the time when the tank become half full.
iii. Find the time when the tank become full.
Problem 1.13. Consider a water thank in the shape of a cylinder of rad ius R = 1 and height H = 1.
There is a drain hole of area A
o
= 0.01m
2
in the bottom of the tank that let water run out of he tank
with the speed Q
o
= 2g
p
A
o
h
p
m
3
/s, where h(t) is the level of the water in the tank and g
=
∼
9.8. If th e
tank is initially full, find the time when the water level drops to half of H.
Problem 1.14. Consider the electrical circuit shown in the following figure.
V
c
CV
S
R
The main tool for analyzing electrical circuit is the Kirchhoff's laws. The mesh law for this circuit
states V
R
+V
c
=V
s
, where V
R
; V
C
are the voltage across the resistance R, and the capacitor C respectively,
and V
s
is the power supply of the circuit. The current-voltage relation across R obeys the law V
R
= Ri,
and the voltage across the capacitor is i = C
dV
c
dt
where i is the electric current in the circuit.
a) Show that the differential equation for V
c
is as follows
RC
dV
c
dt
+ V
c
= V
s
:
b) Verify that the solution to the equation is as follows if V
s
is constant and V
c
(0) = V
0
V
c
(t) = V
s
(1 ¡e
¡t/RC
) + V
0
e
¡t/RC
;
and find lim
t! 1
V
c
(t).
Problem 1.15. Consider the electrical circuit shown in the figure (1.1 ).
R
1
R
2
i
C
i
1
V
c
i
2
Figure 1.1.
a) Show that the differential equation describing V
c
in this case again is
dV
c
dt
+
1
RC
V
c
= 0;
where R =
R
1
R
2
R
1
+ R
2
. (Hint: you need the Kirchhoff's law i = i
1
+ i
2
in this case).
b) Verify that the steady state of V
c
is zero regard less of the initial condition V
c
(0) = V
0
.
Problem 1.16. Consider a water tank that contains V
0
cubic meter of pure water at time t= 0. Suppose
that a current of c
0
gr/lit salty water runs into the tank with speed v
0
lit/s. The content is kept thoroughly
mixed, and simultaneously, the same amount of water runs out of the tank. We would like to determine
the salt concentration in the tank at any instance of time t > 0.
a) Let C(t) denote the total amount of salt in the tank at time t. The conservation princ iple states
dC
dt
= total salt going in¡total s alt going out:
12 Introduction
By the above relation, derive the following differential equation for c(t)
dc
dt
= ¡
v
0
V
0
c +
c
0
v
0
V
0
:
b) Show that the above equation is solved for the following function
c(t) = c
0
(1 ¡e
¡v
0
t/V
0
):
Note that c(0) = 0.
Problem 1.17. Consider a water tank with initial V
0
lit of pure water. Assume that a current of c
0
gr/lit salty water runt into the tank with the speed v
0
lit/sec.
a) Show that the differential equation describing the salt concentration in the tank is
(v
0
t + V
0
)
dc
dt
+ v
0
c = c
0
v
0
:
b) Find the solution to the equation.
Problem 1.18. A mass m = 0.01kg is released from the height h
0
with the initial velocity v
0
=0. Assume
that the air resistance against the mass is f
r
= ¡9.8 ×10
¡6
v
2
.
i. Write down the differential equation for the velocity f unct ion v = v(t) of the m ass (ass ume
g = 9.8m/s
2
).
ii. Find the terminal velocity v¯(the velocity when t !1)
iii. Determine the time t = T when v(T ) = 0.95 v¯.
iv. Determine h
0
such that the velocity of m is exactly equal to 0.95v¯ when it hits the ground.
Problem 1.19. We use the Newton's second law to derive a mathematical model of the motion of a
mass m that is connected to a spring k. Consider the following figure
k
m
x
where x is the displacement of mass m wit h respect to its resting position. The Newton's second
law states
m
d
2
x
dt
2
= f ;
where f is the force applying on the mass. The H o o k 's law states that f = ¡kx, a linear relationship
between the force and t he contraction-stretching of the spring.
a) Derive the differential equation and verify that the obtained equation has two solutions φ
1
=
cos
k
m
q
t
, φ
2
= sin
k
m
q
t
.
b) Now assume that th e ground is not frictions-less and it applies a force to m during its motion on
the surface. The f riction force is f = ¡αv for some constant α. Derive the equation in this case.
1.3 Vector field and the ge ometry of solution
1.3.1 Existence a n d uniqueness
Definition 1.1. Consider the first-order equation
y
0
= f(x; y):
1.3 Vector field and the geometry of solution 13
A continuously differentiable function y = y(x) is called a solution in an open interval I if
y
0
(x) = f (x; y(x)); 8x 2I:
In addition, if an init ial condition y(x
0
) = y
0
is set to the equation, then y(x) is called a
solution to the associated initial value problem if there is an open interval I = (t
0
¡δ; t
0
+ δ)
if y(x
0
) = y
0
and y
0
(x) = f (x; y(x)) for all x 2I.
It is not trivial that an initial value problem possess a solution and that solution is unique.
The following example ex plains the point.
Example 1.2. Consider the following equation
(
y
0
=
y
x
y(0) = 1
:
It is seen that the equation does not possess any solution. In fact, the differential form of
the equation is
dy
y
=
dx
x
;
and thus the general solution is y(x) = cx for arbitrary constant c. Obviously, the obtained
general solution does no t satisfy the initial con dition. Now, consider the following initial
value p roblem
(
y
0
=
y
x
y(0) = 0
:
The equation has infinitely ma ny solution y(x)= cx for arbitrary c. The initial value problem
(
y
0
=
y
x
y(1) = 1
;
has the unique solution y(x) = x.
The existence and uniqueness problem of ordinary differential equations are a deep ques-
tion in the theory of differential equations. We introduce the topic only in an elementary
level in this book.
1.3.2 Vector field of scalar e quations
Consider the first-order equation
dy
dx
= f(y);
for a smooth function y = y(x). Equations that are not explicitly functions of time called
autonomous equations. We can interpret f(y) as a vector field over the y-axis. For example,
consider the following equation
y
0
= 1 ¡y
2
:
14 Introduction
The following figure shows the vector field f(y) = 1 ¡ y
2
over the y-axis
y = 1
y = ¡1
y
Two points y = ±1 where the vector field va nishes are called the equilibrium or cr itical
points of the field. This filed enables us to determine the behavior of the solution qualita-
tively. If y
0
> 1 or ¡1 < y
0
< 1, then we expect
lim
x!1
y(x) = 1:
In this case, the point y = 1 is called an attractor or asymptotically stable equilibrium for
the equation. The point y = ¡1 is called repeller o f the equation.
Definition 1.2. An equilibrium y¯ is called asymptotically stable if there is an interval I =
(y¯ ¡"; y¯+ ") for some " > 0 such that the solution of the following initial value
y
0
= f(y)
y(0) = y
0
2I
; (1.7)
approaches y¯ in long term, that is,
lim
x!1
y(x) = y¯: (1.8)
Theorem 1.2. Let y¯ be an equilibrium of the equation y
0
= f(y) where f is a continuously
differentiable function. The equilibrium is asymptotically stable if f
0
(y¯) <0, and it is unstable
if f
0
(y¯) > 0. The case f
0
(y¯) = 0 is called degenerate case and y¯ can be stable, unstable or
none of them.
See the following figure for a justification of the theorem. If f
0
(y¯) > 0, there is an interval
(y¯ ¡"; y¯ + ") for which y
0
> 0 if y¯+ " > y > y¯ and y
0
< 0 if y¯ ¡" < y < y¯. Therefore, y(x) is
increasing in y¯+" > y > y¯ and de creasing in y¯¡" < y < y¯. Therefore, y¯ is unstable. A similar
argument holds if f
0
(y¯) < 0.
f(y) > 0
f(y) < 0
y: decreasing y
y: increasing
f(y)
f(y) > 0
y: increasing
Figure 1.2.
Problem 1.20. The concavity of the solutions of an a utonomous equation y
0
= f(y) is determined by
d
2
y
dt
2
. D o the concavity a na lysis for the equa tion y
0
= 1 ¡ y
2
. Use the information and draw the graph of
the solution of the equation for y(0) = 0.
1.3 Vector field and the geometry of solution 15
Problem 1.21. By the aid of slope or vector field concept, we can calculate the solutions of differential
equations numerically. Consider the following equation
(
y
0
= 1 ¡ y
2
y(0) = 0
:
We want to find a numerical so lution for y(t) for t 2[0; 1]. Divide the interval into 10 segments. x
0
= 0;
x
1
= 0.1; :::; x
10
= 1.
a) Use y
0
(x
0
) directly by the given equation.
b) Use y
0
(x
0
) by the first order approximation
y(x
1
) = y(x
0
) + y
0
(x
0
)(x
1
¡x
0
):
c) Repeat the above argument and obtain y(x
10
). Use MatLab or other software and draw the true
solution of the equation and the obtained numerical solution.
1.3.3 Vector field of systems
Now consider the following system
8
<
:
dx
dt
= x ¡αx
2
¡0.1xy
dy
dt
= ¡0.5y + 0.1xy
: (1.9)
The vector function F =
x ¡ αx
2
¡ 0.1xy
¡0.5 y + 0.1xy
!
defines a vector field in the (x; y)-plane. This pla ne
is called the phase plane of the system. The following figure shows the vector field for α =0.1.
0 2 4 6 8 10 12
0
1
2
3
4
5
6
7
8
9
10
The e quilibrium point is a point at which
dx
dt
= 0;
dy
dt
= 0. A simple calculation shows that
there are two equilibrium for the system: (0; 0); (5; 5). The solution curves approaches the
second equilibrium in a spiral way. These curves are tangent at everywhere to the vector
field F .
16 Introduction
In general, the solution of a two dimensional autonomous system of the form
8
<
:
dx
dt
= f(x; y)
dy
dt
= g(x; y)
;
is a planar curve γ(t; p
0
) ⊂ R
2
where p
0
= (x
0
; y
0
) is the initial condition of the system
x(0) = x
0
; y(0) = y
0
such that γ(0; p
0
) = p
0
. In addition, the tangent vector γ
0
(t; p
0
) at any
time t is equal to
f(γ(t))
g(γ(t))
.
Problems
Problem 1.22. In the equation (1.9) , assume α = 0 . Use MatLab t o draw the vector field and a few of
trajectories of the system. If you are using MatLab, you can use the streamslice co mmand for trajectories
and quiver for the vector field. The code of the above figure is as follows
[x,y]=meshgrid(0:1:12,0:1:10);
dx=x-0.1*x.^2-0.1*x.*y;dy=-0.5*y+0.1*x.*y;
streamslice(x,y,dx,dy,0.5)
hold on
norm=sqrt(dx.^2+dy.^2);
quiver(x,y,dx./norm,dy./norm,0.5)
axis tight
Problem 1.23. Draw th e vector field and a few of trajectories of the following system
x
0
= y
y
0
= ¡sin(x)
Problem 1.24. For the equation y
0
= x ¡2 y, determine regions in the plane (x; y) where the slope field
is deceasing and increasing respectively. Draw few of the slopes in each region by hand.
Problem 1.25. For a particle that moves under the equation y = x ¡2y, find an approximate path in
the pl ane (x; y) passing through the initial point (0; 0.5). Take h = 0.5 and find point y(0.5), y(1) and
y(1.5). Connect these point and draw the path by hand or use a computer software.
Problem 1.26. Consider a particle moving along the y-axis that its motion obeys the equation y
0
=
1 ¡ y
2
.
a) What is the speed of the particle when y = 3? (note that the velocity is vector and speed is the
length of t his vector).
b) Use the linear approximation formula and find f(h) for h = 0.5 if y(0) = 3. Why t he result is
wrong? Now assume h smaller, say h = 0.2. Find y(0.2), and then use the obtained value to
estimate y(0.4).
c) If y
0
= ¡0.9, find y(0.2) and y(0.4) by the successive linear approximation.
d) Find y(0.2) and y(0.4) if y
0
= ¡1.2.
Problem 1.27. Use a computer software to draw the slope field of the following differential equations
and draw some solution curves. (You can use online applications or use the free software wxMaxima for
it. See the appendix of this chapter for a quick help on wxMaxima)
i. y
0
= sin(y)
ii. y
0
= sin(xy)
iii. y
0
= x ¡e
¡y
1.3 Vector field and the geometry of solution 17
iv. y
0
= e
x¡ y
Problem 1.28. Match the following equations with the given slope fields in the figure (1.3)
i. y
0
= ¡y
ii. y
0
= x + y
iii. y
0
= x ¡y
2
iv. y
0
= 1 ¡ y
2
-2 -1 1 2
-2
-1
1
2
-2 -1 1 2
-2
-1
1
2
-2 -1 1 2
-2
-1
1
2
-2 -1 0 1 2
-2
-1
0
1
2
Figure 1.3.
Problem 1.29. Consider the equation y
0
= y ¡y
3
.
i. Find equilibrium points of the equation.
ii. Do stability analysis for obtained equilibrium points.
iii. Draw the integral curve passing through
¡
0;
1
4
.
iv. Verify that the function φ(x) =
1
1 + 15 e
¡2x
p
is a solution to the equation satisfying th e condition
y(0)=
1
4
. Use a computer software to draw this function and compare it with the result of part (iii).
Problem 1.30. Consider the equation y
0
= y
2
¡y
3
.
i. If y(0) =
1
2
find lim
x!1
y(t) an d draw th e solution,
ii. If y(0) =
¡1
2
find lim
x!1
y(t) and draw the solution,
iii. If y( 0) = 2 find lim
x!1
y(t) an d draw the solution.
Problem 1.31. For each of the following equation s, find equilibrium points, and do stability analysis
and draw the solution curve
i. y
0
= ¡y
2
+ 3y ¡2, y(0) = 3/2
18 Introduction
ii. y
0
= sin (1 ¡ y), y(0) = π /2
iii. y
0
= y(1 ¡y)(2 ¡ y), y(0) = 3
iv. y
0
= e
y
y(1 ¡ y), y(0) = 1/2
v. y
0
= (1 ¡ y
p
) y, y(0) = 1/2
Problem 1.32. Consider the equation
y
0
+ xy = x
2
+ 1:
i. At what angel does the graph of the sol ution cut the y-axis?
ii. If the graph of the solution intersects t-axis, show that the intersection angle θ lies in the range
π
4
;
π
2
.
Problem 1.33. Consider the problem
(
(1 + y
2
)y
0
= 1 + x
2
y(0) = ¡1
i. Find the slope of the solution curve y = φ(x) at x = 0.
ii. Prove that the solution curve y = φ (x ) crosses the x-axis.
iii. Prove th at the solution curve y = φ(x) crosses the tx-axis with the angle greater than
π
4
.
Problem 1.34. Consider the initial value problem
y
0
= sin(x)sin(y) + 1
y(0) = ¡1
:
i. Show that the solution curve intersects the x-axis.
ii. Show that the solution curve intersects the x-axis only at one point.
iii. Find the angle of the intersection.
1.4 Higher o rder equations
1.4.1 The general form
The general form o f a second-order differential equation for a function y = y(x) is
F (x; y; y
0
; y
00
) = 0; (1.10)
where y
00
=
d
2
y
dx
2
. An initial value problem for a second-order equation is of the following form
8
<
:
F (x; y; y
0
; y
00
) = 0
y(x
0
) = y
0
y
0
(x
0
) = y
1
:
A differential equation of order n is
F (x; y; y
0
; :::; y
(n)
) = 0; (1.11)
1.4 Higher order equations 19
where y
(n)
=
d
n
y
dx
n
. An initial value problem of order n is of the following form
8
>
>
>
>
>
>
>
>
>
>
<
>
>
>
>
>
>
>
>
>
>
:
F (x; y; y
0
; :::; y
(n)
) = 0
y(x
0
) = y
0
y
0
(x
0
) = y
1
·
·
·
y
(n¡ 1)
(x
0
) = y
n¡1
:
If function F is linear with respect to y and its derivatives, then the differential equation is
called linear. The general form of a linear first order eq uation is
p(x)y
0
+ q(x)y = r(x):
A second order lin ear equation has the following form
p(x)y
00
+ q(x)y
0
+ z(x) y = r(x):
If r(x) is identically zero, the equation is ca lle d linear homogeneous equation. A linear
differential equation of order n is of the following general form
a
n
(x) y
(n)
+ a
n¡1
y
(n¡ 1)
+ ···+ a
0
(x)y = r(x):
1.4.2 Nonlinear equa tions
A differential equation that is not linear is called nonlinear different ial equation. Nonlinear
equation are extremely difficult to solve and sometime show strange or chaotic behavior.Con-
sider the pendulum shown in the f o llowing figure
l θ l
m
m
The motion of m fo llows the following equation
d
2
θ
dt
2
+
g
l
sin(θ) = 0:
The equation is non-linear due to the no n-linear term sin(θ). But for θ small, si n(θ) ≈θ, and
physicists tend to rewrite the equation as the following linear one
d
2
θ
dt
2
+
g
l
θ = 0:
The general solution of the above linear equation is
θ(t) = c
1
cos
g
l
r
t
!
+ c
2
sin
g
l
r
t
!
;
20 Introduction
while the solution to the original equation is complicated. Let us continue the example and
apply an external force on the mass of the form f = sin(2t). The original equation with this
external force reads
d
2
θ
dt
2
+
g
l
sin(θ) = sin(2t):
The following figure shows the so lution of both equations in the (θ; θ
0
)-plane (this is called
the phase plane). The left figure is for the nonlinear one, while the right figure is f o r the
linear one. As it is observed, the motion of the nonlinear forc ed pendulum is chaotic, while
the solution of linear equation is very regular.
θ
θ
′
θ
θ
′
Problems
Problem 1.35. For each of the following equations, determine the order of equa tion, if it is linear or
non-linear and also if it is autonomous or non-autonomous:
i. y
0
+ xy = e
y
ii. y
00
+ yy
0
= 0
iii. y
(4)
+ y (y
000
)
5
= sin(t)
iv. y
00
+ sin(y) = 0
v. x
2
y
00
+ 2xy
0
+ 2y = e
x
Problem 1.36. Consider the linear equation y
0
¡tan(x)y = 2x sec(x).
a) Verify that the function y = x
2
sec(x) is a solution to the equation.
b) Now consider the initial condition y(0) = 1. Obviously, the solution y = x
2
sec(x) does not satisfies
the initial condition. Verify that the function y
h
= c sec(x) solves the equation y
0
¡tan(x) y = 0
for arbitrary constant c, and thus conclude that the function y = c sec(x) + x
2
sec(x) satisfies also
the original equation.
c) Now, find a solution y that satisfies the equation and the initial condition y(0) = 1.
Problem 1.37. Consider the linear equation y
0
+ 2xy = x. Obviously y =
1
2
is a solution to the equation,
however, the ha s other solutions too.
a) Show that the function y
h
= ce
x
2
solves the homogeneous equation y
0
+ 2xy = 0 and conclude that
the function y = ce
x
2
+
1
2
solves the original equation.
b) Find and equation that satisfies y(0) = 0.
1.4 Higher order equations 21
Problem 1.38. Consider the linear equation y
0
+ p( x)y = 0.
a) Verify that if y is a solution to the linear equation then y
h
= cy(x) is also the sol ution for all
constants c.
b) This property does not hold for non-linear equations in general. Consider the equation y
0
+ y
2
= 0.
Verify that y =
1
x
is a solution to the equation. Is the function y =
c
x
the solution to the equation
for all c?
Problem 1.39. Consider the linear equation y
00
¡2y
0
+ y = 0.
a) Verify that functions y
1
= e
x
; y
2
= xe
x
are solutions to the equation.
b) Verify tha t for any constant c
1
; c
2
, the function y = c
1
y
1
+ c
2
y
2
is a solution to the equation as well.
Problem 1.40. Verify that the function y = c
1
x + c
2
x
2
is a solution to the equation
x
2
y
00
¡2xy
0
+ 2y = 0:
Find the so lution to the equation satisfying initial conditions y(1) = 1; y
0
(1) = 0. Is it possible to find a
solution to the equation satisfying the condition y(0) = 1, y
0
(0) = 0?
Problem 1.41. Verify that the function y =c
1
cosh(x)+ c
2
sinh(x) is a solution to the equation y
00
¡y =0.
Find the solution to the equation satisfying initial conditions y(0) = y
0
(0) = 1.
Problem 1.42. Consider the equation y
0
= y ¡y
3
:
i. Verify tha t the function y
2
¡
1
1 + ce
¡2x
= 0 for arbitrary constant c is a solution to the equation.
ii. Find an explicit solution y = φ(x) for the equation satisfying the initial condition y(0) = 1. Find
lim
x!1
y(x) and use a computer software to draw the solution in its domain of definition.
iii. Find an explicit solution y = φ(x) to the equation satisfying the initial condition y(0) = ¡1. Find
lim
x!1
y(x) and draw the solution in its domain of definition.
iv. Find an explicit solution to the equation satisfying the initial condition y(0) = 2 and find the
domain of the definition for the solution.
v. Find the solution to the equation satisfying the initial condition y(0) = 1.
Problem 1.43. Consider the equation (1 ¡x
2
)y
0
¡2xy = 0.
a) Verify that the function y =
c
1 ¡ x
2
is a solution to the equation f or arbitrary constant c.
b) Determine c if y(0) = 1 and obtain the maximum interval that this solution extends.
c) Repeat part b) if y(2) = 1.
Problem 1.44. Here we see a simple initial value problem without any possible solution. Consider the
equation (x + 1)y
0
= 1 + y. Verify that the function y = c(x + 1) ¡1 solves the equation for all constant
c, and conclude that there is no solution to th e initial value problem (1 + x)y
0
= 1 + y, y(¡1) = 1.
Problem 1.45. Consider the equation y
0
= cos(y).
a) Verify that the implicit function sec(y) +t an (y) = ce
x
solves the equation for arbitrary constant c.
b) Determine c if y(0) = 0 and use the implicit function theorem to show that there is an explicit
solution for the equation. Can you determine the possib le interval the solution extend?
Problem 1.46. Find a solution to the problem
(
x cos(xy) y
0
+ cos(xy) y = e
x
y(1) =
π
2
;
and determine the domain of definition for the solution. Draw the solution in its domain. (H int: the left
hand side is the derivative of the function sin(xy)).
22 Introduction
Problem 1.47. Try to find a solution to the following equation
(cos(y) ¡ y sin(y)) y
0
= 2x
y(0) = y
0
:
Draw integral curves for y
0
= ¡2; ¡1; ¡0.5; 0.5 in the interval ¡3 < x < 3, ¡6 ≤y ≤4. Find the solution
that passes through the point (0; 1).
Problem 1.48. Verify that the following functions solve the heat equation defined above.
u = e
¡D!t
cos( !
p
x); u = e
¡D!t
sin( !
p
x):
Problem 1.49. The partial differential equation
@
2
u
@x
2
+
@
2
u
@y
2
= 0;
is called the harmonic equation for a function u = u(x; y). Verify that the following functions are solution
to the harmonic equation
u = e
x
cos(y); y = cosh(x) sin(y); y = x
2
¡y
2
:
1.4 Higher order equations 23
Appendix A
Basic commands in wxMaxima
We review some basic commands in wxMaxima math software. The complete explanation
can be found in the tutorial of WxMaxima.
To integrate a function f (x)
integrate(f,x);
For definite integral in the range [a; b]
integrate(f,x,a,b);
The following simple c od e plots the function f from x = x
0
to x = x
1
:
plot2d(f,[x,x0,x1]);
Let us integrate the function f (x) = x s in(nx) in the range x 2[0; π] fo r n nonzero integers
declare(n,integer);
integrate(x*sin(n*x),x,0,%pi]);
The summation sum:
sum(n,n,1,10);
adds numbers from 1 to 10. Similarly the command
sum(2*(1-cos(n*%pi))*sin(n*%pi*x)/(n*%pi),n,1,5)
adds the given function from 1 to 5.
For the plot of a function f (x) in the range [a; b]
plot2d(f,[x,a,b]);
For example for f = sin(x) from x = ¡π to π we can write
plot2d(sin(x),[x,%pi,%pi]);
You can save the o utput in a p s, pdf or png as foll ows
plot2d(f,[x,x0,x1],[ps_file,"filename.ps"]);
plot2d(f,[x,x0,x1],[pdf_file,"filename.pdf"]);
plot2d(f,[x,x0,x1],[png_file,"filename.png"]);
For multiple functions f
1
; f
2
; ··· we can draw them in a same coordinate
plot2d([f1,f2,...],[x,x0,x1]);
For example to draw f
1
= e
0.5x
, f
2
= log(jxj) a nd f
3
= 0.5e
x
cos(2x) in ¡1 ≤x ≤4, write
plot2d([%e^(0.5*x),log(abs(x)),0.5*%e^x*cos(2*x)],[x,-1,4]);
If you want to hide the le g end, just use the command
plot2d([%e^(0.5*x),log(abs(x)),0.5*%e^x*cos(2*x)],[x,-1,4],
[legend,false]);
To draw a parametric curve (x(t); y(t)), from t = t
0
to t = t
1
use
plot2d([parametric,x(t),y(t),[t,t0,t1]]);
For example, the circle (cos(t); sin(t)) for t = 0 to t = 2π can be plotted by
plot2d([parametric,cos(t),sin(t),[t,0,2*%pi]]);
25
To make the axes ratio same, use the command [yx_ratio,1]. For example
plot2d([parametric,cos(t),sin(t),[t,0,2*%pi]],[yx_ratio,1]);
To draw a list data us the following command
plot2d([discrete, data]);
For example to plot the list [ [¡3; 9 ]; [¡ 2 ; 4]; [¡1; 1]; [0; 0]; [1; 1]; [2; 4]; [3; 9]] use
data:makellist([i,i^2],i,-3,3);
plot2d([discrete, data]);
To draw 3D plots
plot3d(f,[x,x0,x1],[y,y0,y1]);
For example
plot3d(%e^(x+y)*sin(x*y),[x,0,%pi],[y,0,%pi]);
To draw an implicit function f(x; y) = 0, load first the package as
load(implicit_plot);
and then draw by
implicit_plot(f(x,y)=0,[x,x0,x1],[y,y0,y1]);
For example the following command draw a circle with radius 2:
implicit_plot(x^2+y^2=4,[x,-2,2],[y,-2,2]);
To draw the slope field (or vector field) of a differential equation, you should fir st load
the package plotdf as
load(plotdf);
and then to draw the slope field. For example to draw the field of the equation y
0
= sin(y)
just write
plotdf(sin(y),[x,x0,x1],[y,y0,y1]);
When you click on the plot page, a trajectory is appeared passing through the click point.
In order to draw the vector field [f(x; y); g(x; y)] just write
plotdf([f,g],[x,x0,x1],[y,y0,y1]);
For example the vector field V = (¡y; x) in ¡1 ≤x ≤1 and ¡1 ≤y ≤1 can be plotted as
plotdf([-y,x],[x,-1,1],[y,-1,1]);
To draw a trajectory starting at, say (0.5; 0.5) write
plotdf([-y,x],[x,-1,1],[y,-1,1],[trajectory_at,0.5,0.5]);
To sol ve an initial val ue problem numerically, use the Runge-Kutta method. The
command
rk(f(x,y),y,y0,[x,x0,x1,xstep]);
solves the i.v.p. y
0
= f(x; y), y(x
0
) = y
0
for x
0
≤x ≤x
1
for steps xstep.
For example, the i.v.p. y
0
= cos(2x) + ye
¡y
, y(0) = 1 is solved numerically for 0 ≤x ≤5 as
rk(cos(2*x)+y*%e^(-y),y,1,[x,0,5,0.1]);
To draw the integral curve o f the above i.v.p. use the following command
data:rk(cos(2*x)+y*%e^(-y),y,1,[x,0,5,0.1]);
plot2d([discrete,data]);
For a system of equations y
0
= f(x; y; z), z
0
= g(x; y; z) do as follows.
data:rk([f,g],[y,z],[y0,z0],[x,x0,x1,xstep]);
plot2d([discrete,makelist([p[1],p[2]],p,data)]);
that above command plot y(x). To draw z(x) use
plot2d([discrete,makelist([p[1],p[3]],p,data)]);
and to draw the phase portrait (y(x); z(x)) use
26 Basic commands in wxMaxima
plot2d([discrete,makelist([p[2],p[3]],p,data)]);
For constructing an approximate integral curve for a given vector field V (¡y; x), you can
do as follows. Fist define the vector field
V(x,y):=(-y,x);
Then initialize the step size, say for h = 0.02
h:0.2;
and the initial point, say for p
0
= (1; 0)
p[0]:[1,0];
Define the recursive formula
p[n]:=p[n-1]+h*V(p[n-1][1],p[n-1][2]);
and then make a list, for example for 200 entries
data:makelist(p[n],n,1,200);
Now you can draw it as
plot2d([discrete,data]);
Let us solve the following i.v.p.
y
00
+ cos(x) y
0
+ xy = e
x
y(0) = 1; y
0
(0) = ¡1
:
Step 1. First we transform the equation into a system of two first order equations a s
follows. If we take y = y
1
then the following system is equivalent to the equation:
y
1
0
= y
2
y
2
0
= e
x
¡xy
1
¡cos(x) y
2
:
Step 2. The initial conditions:
y
1
(0) = 1; y
2
(0) = ¡1
Step 3. Use the Runge-Kutta built in method in WxMaxima. The general form is
rk([system],[state variable],[initial conditions],[x,0,x1,step-size])
For our system and x
1
= 2 with the step size 0.1 it will be
rk([y2,-exp(x)-x*y1-cos(x)*y2],[y1,y2],[1,-1],[x,0,2,0.1])
Note that for smaller step size the solution is more accurate.
The output is a list of three entries as
[[0; y
1
(0); y
2
(0)]; [0.1; y
1
(0.1); y
2
(0.1)]; [0.2; y
1
(0.2); y
2
(0.2)]; :::]
In order to plot th e output list [:::; [x; y
1
(x)]; :::] we do as follows. First, let us assign a handle
to the solution
sol: rk([y2,-exp(x)-x*y1-cos(x)*y2],[y1,y2],[1,-1],[x,0,2,0.1]);
Then make a list as follows
list:makelist([p[1],p[2]],p,sol);
Finally plot by the command plot2d in WxMaxima
plot2d([discrete,list]);
If you want plot the output an a function f (x) is a same coordinate you can do as
plot2d([f(x),[discrete,list]],[x,0,2]);
Basic commands in wxMaxima 27
