Chapter 3: 2D linear systems (including the harmonic oscillator)

>    diff(x(t),t) = a(t)*x(t) + b(t)*y(t) + f(t);
diff(y(t),t) = c(t)*x(t) + d(t)*y(t) + g(t);

diff(x(t),t) = a(t)*x(t)+b(t)*y(t)+f(t)

diff(y(t),t) = c(t)*x(t)+d(t)*y(t)+g(t)

Constant coefficient systems

If a, b, c, and d are all constants, the system has constant coefficients.

Two-dimensional linear constant coefficient systems are of the form

>    diff(x(t),t) = a*x(t) + b*y(t) + f(t);
diff(y(t),t) = c*x(t) + d*y(t) + g(t);

diff(x(t),t) = a*x(t)+b*y(t)+f(t)

diff(y(t),t) = c*x(t)+d*y(t)+g(t)

Homogeneous systems

If f and g are both zero, the system is homogeneous .
Two-dimension linear homogeneous systems are of the form

>    diff(x(t),t) = a(t)*x(t) + b(t)*y(t);
diff(y(t),t) = c(t)*x(t) + d(t)*y(t);

diff(x(t),t) = a(t)*x(t)+b(t)*y(t)

diff(y(t),t) = c(t)*x(t)+d(t)*y(t)

Constant coefficient homogeneous systems part 1:
   componentwise derivation of the characteristic equation

>    DE1 := diff(x(t),t) = a*x(t)+b*y(t): DE1;
DE2 := diff(y(t),t) = c*x(t)+d*y(t): DE2;

diff(x(t),t) = a*x(t)+b*y(t)

diff(y(t),t) = c*x(t)+d*y(t)

>    xguess := x0*exp(lambda*t): x(t) = xguess;
yguess := y0*exp(lambda*t): y(t) = yguess;

x(t) = x0*exp(lambda*t)

y(t) = y0*exp(lambda*t)

>   

>    E1 := simplify(subs(x(t)=xguess,y(t)=yguess,DE1)): E1;
E2 := simplify(subs(x(t)=xguess,y(t)=yguess,DE2)): E2;

x0*lambda*exp(lambda*t) = exp(lambda*t)*(a*x0+b*y0)

y0*lambda*exp(lambda*t) = exp(lambda*t)*(c*x0+d*y0)

>    E3 := simplify(E1/exp(lambda*t)): E3;
E4 := simplify(E2/exp(lambda*t)): E4;

x0*lambda = a*x0+b*y0

y0*lambda = c*x0+d*y0

>    E5 := (a-lambda)*x0 + b*y0 = 0;
E6 := c*x0 + (d-lambda)*y0 = 0;

E5 := (a-lambda)*x0+b*y0 = 0

E6 := c*x0+(d-lambda)*y0 = 0

>    E7 := simplify(E5*(d-lambda)-E6*b);

E7 := d*a*x0-d*x0*lambda-lambda*a*x0+x0*lambda^2-b*c*x0 = 0

>    E7 := factor(E7);

E7 := x0*(d*a-d*lambda-lambda*a+lambda^2-c*b) = 0

>    coeff(lhs(E7),x0) <> 0;

d*a-d*lambda-lambda*a+lambda^2-c*b <> 0

>    E8 := factor(simplify(E6*(a-lambda)-E5*c));

E8 := y0*(d*a-d*lambda-lambda*a+lambda^2-c*b) = 0

>    coeff(lhs(E8),y0) <> 0;

d*a-d*lambda-lambda*a+lambda^2-c*b <> 0

>    lambda^2 - (a+d)*lambda + a*d - c*b = 0;

lambda^2-(a+d)*lambda+d*a-c*b = 0

>   

Constant coefficient homogeneous systems part 2:
    the eigenproblem

>    with(linalg):
colvector:=proc(rowvec) matrix(nops(rowvec),1,rowvec) end;
u := colvector([x(t),y(t)]);
dudt := map(f->diff(f,t),u);
 

Warning, the name adjoint has been redefined



Warning, the protected names norm and trace have been redefined and unprotected

colvector := proc (rowvec) matrix(nops(rowvec),1,rowvec) end proc

u := matrix([[x(t)], [y(t)]])

dudt := matrix([[diff(x(t),t)], [diff(y(t),t)]])

>    DE1 := diff(x(t),t) = a*x(t)+b*y(t): DE1;
DE2 := diff(y(t),t) = c*x(t)+d*y(t): DE2;

diff(x(t),t) = a*x(t)+b*y(t)

diff(y(t),t) = c*x(t)+d*y(t)

>    DE := eval(dudt) = colvector([a*x(t)+b*y(t),c*x(t)+d*y(t)]): DE;

matrix([[diff(x(t),t)], [diff(y(t),t)]]) = matrix([[a*x(t)+b*y(t)], [c*x(t)+d*y(t)]])

>    A:=matrix([[a,b],[c,d]]);

A := matrix([[a, b], [c, d]])

>    eval(A)&*eval(u) = evalm(A&*u);

`&*`(matrix([[a, b], [c, d]]),matrix([[x(t)], [y(t)]])) = matrix([[a*x(t)+b*y(t)], [c*x(t)+d*y(t)]])

>    DE := dudt = A&*u: DE;
DE := eval(dudt) = eval(A)&*eval(u): DE;

dudt = `&*`(A,u)

matrix([[diff(x(t),t)], [diff(y(t),t)]]) = `&*`(matrix([[a, b], [c, d]]),matrix([[x(t)], [y(t)]]))

>    u = colvector([x0*exp(lambda*t),y0*exp(lambda*t)]);

u = matrix([[x0*exp(lambda*t)], [y0*exp(lambda*t)]])

>    DE0 := map(g->simplify(subs(x(t)=x0*exp(lambda*t),y(t)=y0*exp(lambda*t),g)),DE): DE0;

matrix([[diff(x0*exp(lambda*t),t)], [diff(y0*exp(lambda*t),t)]]) = `&*`(matrix([[a, b], [c, d]]),matrix([[x0*exp(lambda*t)], [y0*exp(lambda*t)]]))

>    u0:= colvector([x0,y0]):
eval(A)&*eval(u0) = lambda*eval(u0);
A&*u0 = lambda*u0;

`&*`(matrix([[a, b], [c, d]]),matrix([[x0], [y0]])) = lambda*matrix([[x0], [y0]])

`&*`(A,u0) = lambda*u0

>    I2 := matrix([[1,0],[0,1]]): I2 = evalm(I2);

I2 = matrix([[1, 0], [0, 1]])

>    (A-lambda*I2)&*u0 = 0;
evalm(A-lambda*I2)&*eval(u0) = 0;

`&*`(A-lambda*I2,u0) = 0

`&*`(matrix([[a-lambda, b], [c, d-lambda]]),matrix([[x0], [y0]])) = 0

>    (a-lambda)*(d-lambda) - c*b = 0;

(a-lambda)*(d-lambda)-c*b = 0

>    chareqn := lambda^2 - (a+d)*lambda + a*d - b*c = 0: chareqn;

lambda^2-(a+d)*lambda+d*a-c*b = 0

>   

Eigenvalues and eigenvectors

>    A := matrix([[a,b],[c,d]]):
I2 := matrix([[1,0],[0,1]]):
colvector := proc(expr) matrix(nops(expr),1,expr) end:
u0 := colvector([x0,y0]):

>    A&*u0 = lambda*u0;
eval(A)&*eval(u0) = lambda*eval(u0);
evalm(A-lambda*I2)&*eval(u0) = 0;

`&*`(A,u0) = lambda*u0

`&*`(matrix([[a, b], [c, d]]),matrix([[x0], [y0]])) = lambda*matrix([[x0], [y0]])

`&*`(matrix([[a-lambda, b], [c, d-lambda]]),matrix([[x0], [y0]])) = 0

>    (a-lambda)*(d-lambda)-c*b = 0;

(a-lambda)*(d-lambda)-c*b = 0

>    evalm(A-lambda[1]*I2)&*eval(u0) = 0;

`&*`(matrix([[a-lambda[1], b], [c, d-lambda[1]]]),matrix([[x0], [y0]])) = 0

>    evalm(A-lambda[2]*I2)&*eval(u0) = 0;

`&*`(matrix([[a-lambda[2], b], [c, d-lambda[2]]]),matrix([[x0], [y0]])) = 0

for a new pair of values x0 and y0. Let v2 be the column vector formed from these x0, y0.

>    A := matrix([[1,4],[4,1]]);

A := matrix([[1, 4], [4, 1]])

>    eval(A)&*eval(u0) = lambda*eval(u0);
evalm(A-lambda*I2)&*eval(u0) = 0;

`&*`(matrix([[1, 4], [4, 1]]),matrix([[x0], [y0]])) = lambda*matrix([[x0], [y0]])

`&*`(matrix([[1-lambda, 4], [4, 1-lambda]]),matrix([[x0], [y0]])) = 0

The characteristic equation is then

>    chareqn :=
(A[1,1]-lambda)*(A[2,2]-lambda)-A[2,1]*A[1,2] = 0;

chareqn := (1-lambda)^2-16 = 0

>    lam := solve(chareqn,lambda);

lam := 5, -3

>    lambda[1]:='lambda[1]':
evalm(A-lambda[1]*I2)&*eval(u0) = 0;
evalm(A-lam[1]*I2)&*eval(u0) = 0;
evalm((A-lam[1]*I2)&*u0 = 0);

`&*`(matrix([[1-lambda[1], 4], [4, 1-lambda[1]]]),matrix([[x0], [y0]])) = 0

`&*`(matrix([[-4, 4], [4, -4]]),matrix([[x0], [y0]])) = 0

matrix([[-4*x0+4*y0], [4*x0-4*y0]]) = 0

>    v1 := colvector([1,-1]);

v1 := matrix([[1], [-1]])

is an eigenvector corresponding to the eigenvalue

>    lambda[1] = lam[1];

lambda[1] = 5

>    evalm(A&*v1);  evalm(lam[1]*v1);

matrix([[-3], [3]])

matrix([[5], [-5]])

These agree, as they should..

>    lambda[2]:='lambda[2]':
evalm(A-lambda[2]*I2)&*eval(u0) = 0;
evalm(A - lam[2]  *I2)&*eval(u0) = 0;
evalm((A - lam[2]  *I2)&*u0 = 0);

`&*`(matrix([[1-lambda[2], 4], [4, 1-lambda[2]]]),matrix([[x0], [y0]])) = 0

`&*`(matrix([[4, 4], [4, 4]]),matrix([[x0], [y0]])) = 0

matrix([[4*x0+4*y0], [4*x0+4*y0]]) = 0

>    v2 := colvector([1,1]);

v2 := matrix([[1], [1]])

>    lambda[2] = lam[2];

lambda[2] = -3

>    evalm(A&*v2); evalm(lam[2]*v2);

matrix([[5], [5]])

matrix([[-3], [-3]])

These agree, as they should.

Solution of a constant coefficient homogenous system with
  coefficient matrix having real unequal eigenvalues

>    A := matrix([[1,4],[4,1]]):
DE := colvector(diff([x(t),y(t)],t)) = evalm(A)&*colvector([x(t),y(t)]): DE;

matrix([[diff(x(t),t)], [diff(y(t),t)]]) = `&*`(matrix([[1, 4], [4, 1]]),matrix([[x(t)], [y(t)]]))

>    lam := 'lam': lam[1] := -3: lam[2] := 5:
lambda[1] := 'lambda[1]': lambda[2]:='lambda[2]':
lambda[1] = lam[1]; v1 := colvector([1,-1]);
lambda[2] = lam[2]; v2 := colvector([1,1]);

lambda[1] = -3

v1 := matrix([[1], [-1]])

lambda[2] = 5

v2 := matrix([[1], [1]])

>    u1 := v1*exp(lambda[1]*t);
u1:= eval(v1)*exp(lam[1]*t);
u2 := v2*exp(lambda[2]*t);
u2 := eval(v2)*exp(lam[2]*t);

u1 := v1*exp(lambda[1]*t)

u1 := matrix([[1], [-1]])*exp(-3*t)

u2 := v2*exp(lambda[2]*t)

u2 := matrix([[1], [1]])*exp(5*t)

>    u_gen := c1*u1 + c2*u2;
u_gen := c1*evalm(u1) + c2*evalm(u2);
u_gen := evalm(c1*u1 + c2*u2);

u_gen := c1*matrix([[1], [-1]])*exp(-3*t)+c2*matrix([[1], [1]])*exp(5*t)

u_gen := c1*matrix([[exp(-3*t)], [-exp(-3*t)]])+c2*matrix([[exp(5*t)], [exp(5*t)]])

u_gen := matrix([[c1*exp(-3*t)+c2*exp(5*t)], [-c1*exp(-3*t)+c2*exp(5*t)]])

>   

>    u1_x :=evalm(u1)[1,1]: u2_x := evalm(u2)[1,1]:
u1_y :=evalm(u1)[2,1]: u2_y := evalm(u2)[2,1]:
plot([u1_x,u2_x],t=0..0.5,title=`x-component of u1 in red, of u2 in green`);
plot([u1_y,u2_y],t=0..0.5,title=`y-component of u1 in red, of u2 in green`);

[Maple Plot]

[Maple Plot]

>    u1xyplot := plot([u1_x,u1_y,t=0..0.5],color=red):
u2xyplot := plot([u2_x,u2_y,t=0..0.5],color=green):
with(plots): display({u1xyplot,u2xyplot},labels=[`x`,`y`],title=`phase plane plots: u1(red), u2(green)`);

[Maple Plot]

>   

Computing eigenvalues and eigenvectors with Maple

>    with(linalg):

>    A := matrix([[1,4],[4,1]]);
charpoly(A,lambda);
eigenvals(A);
eigenvects(A);

A := matrix([[1, 4], [4, 1]])

lambda^2-2*lambda-15

5, -3

[5, 1, {vector([1, 1])}], [-3, 1, {vector([-1, 1])}]

>   

>    colvector([-1,1]);

matrix([[-1], [1]])

>    eval(v1);

matrix([[1], [-1]])

>   

Here's an exercise for you:
   Compute the eigenvalues and eigenvectors of the 2 by 2 identity
matrix I2 above, by using eigenvects(I2).

>   

Solutions of a linear constant coefficient homogenous system
with complex eigenvalues: the harmonic oscillator

>    A := matrix([[0,1],[-1,0]]):
colvector := proc(expr) matrix(nops(expr),1,expr) end:
u := colvector([x(t),y(t)]):
dudt := colvector([diff(x(t),t),diff(y(t),t)]):
dudt = A&*u;
eval(dudt) = eval(A)&*eval(u);

dudt = `&*`(A,u)

matrix([[diff(x(t),t)], [diff(y(t),t)]]) = `&*`(matrix([[0, 1], [-1, 0]]),matrix([[x(t)], [y(t)]]))

>    uguess := eval(u0)*exp(lambda*t): eval(u) = uguess;

matrix([[x(t)], [y(t)]]) = matrix([[x0], [y0]])*exp(lambda*t)

>    evalm(matrix([[A[1,1]-lambda,  A[1,2]],
              [A[2,1], A[2,2]-lambda]]))&*eval(u0) = 0;

`&*`(matrix([[-lambda, 1], [-1, -lambda]]),matrix([[x0], [y0]])) = 0

>    chareqn :=
(A[1,1]-lambda)*(A[2,2]-lambda)-A[1,2]*A[2,1] = 0;
lam := solve(chareqn,lambda);

chareqn := lambda^2+1 = 0

lam := I, -I

>   


To find an eigenvector [x0,y0] corresponding to the first eigenvalue ,
we consider the eigenproblem with lambda = I   :

>    lambda[1]:='lambda[1]': x0:='x0': y0:='y0':
evalm(A-lambda[1]*I2)&*eval(u0) = 0;
evalm(A-   I    *I2)&*eval(u0) = 0;
evalm((A-   I    *I2)&*u0) = 0;
eqn1 := -I*x0 + y0 = 0;
eqn2 := -x0 - I*y0 = 0;

`&*`(matrix([[-lambda[1], 1], [-1, -lambda[1]]]),matrix([[x0], [y0]])) = 0

`&*`(matrix([[-I, 1], [-1, -I]]),matrix([[x0], [y0]])) = 0

matrix([[-I*x0+y0], [-x0-I*y0]]) = 0

eqn1 := -I*x0+y0 = 0

eqn2 := -x0-I*y0 = 0

>    x0 := 1; y0 := solve(eqn1,y0);
v1 := colvector([x0,y0]);

x0 := 1

y0 := I

v1 := matrix([[1], [I]])

>    eval(eqn2);

0 = 0

>    t:='t':
u1 := v1*exp(lambda[1]*t);
u1 := evalm(v1)*exp(I*t);

u1 := v1*exp(lambda[1]*t)

u1 := matrix([[1], [I]])*exp(t*I)

>    eIt := cos(t)+I*sin(t):
u1 := evalm(v1)*(eIt);

u1 := matrix([[1], [I]])*(cos(t)+sin(t)*I)

>    u1 := evalm(v1*eIt);

u1 := matrix([[cos(t)+sin(t)*I], [(cos(t)+sin(t)*I)*I]])

>    u1 := map(f->combine(f),u1);

u1 := matrix([[cos(t)+sin(t)*I], [cos(t)*I-sin(t)]])

>    w1 := colvector([cos(t),-sin(t)]);
w2 := colvector([sin(t),cos(t)]);

w1 := matrix([[cos(t)], [-sin(t)]])

w2 := matrix([[sin(t)], [cos(t)]])

>    map(f->diff(f,t),w1) = evalm(A&*w1);
map(f->diff(f,t),w2) = evalm(A&*w2);

matrix([[-sin(t)], [-cos(t)]]) = matrix([[-sin(t)], [-cos(t)]])

matrix([[cos(t)], [-sin(t)]]) = matrix([[cos(t)], [-sin(t)]])

>    u_gen = c1*w1 + c2*w2;
u_gen = evalm(c1*w1 + c2*w2);

u_gen = c1*w1+c2*w2

u_gen = matrix([[c1*cos(t)+c2*sin(t)], [-c1*sin(t)+c2*cos(t)]])

where c1 and c2 are arbitrary real constants.

>   

>    w1_x := w1[1,1]:  w2_x := w2[1,1]:
w1_y := w1[2,1]:  w2_y := w2[2,1]:
plot([w1_x,w2_x],t=0..10,labels=[`t`,`x`],title=`x-component of w1 in red, of w2 in green`);
plot([w1_y,w2_y],t=0..10,labels=[`t`,`y`],title=`y-component of w1 in red, of w2 in green`);

[Maple Plot]

[Maple Plot]

>    w1xyplot := plot([w1_x,w1_y,t=0..1],-2.7..2.7,-2..2,color=red):
w2xyplot := plot([w2_x,w2_y,t=0..1],color=green):
with(plots): display({w1xyplot,w2xyplot},labels=[`x`,`y`],title=`phase plane plots: w1(red), w2(green)`);

[Maple Plot]

>   

Solutions of a linear constant coefficient homogenous system
with complex eigenvalues, in which the solutions decay exponentially

>    A := matrix([[-1,4],[-1,-1]]):
colvector := proc(expr) matrix(nops(expr),1,expr) end:
u := colvector([x(t),y(t)]):
dudt := colvector([diff(x(t),t),diff(y(t),t)]):
dudt = A&*u;
evalm(dudt) = evalm(A)&*evalm(u);

dudt = `&*`(A,u)

matrix([[diff(x(t),t)], [diff(y(t),t)]]) = `&*`(matrix([[-1, 4], [-1, -1]]),matrix([[x(t)], [y(t)]]))

>    uguess := evalm(u0)*exp(lambda*t): evalm(u) = uguess;

matrix([[x(t)], [y(t)]]) = matrix([[x0], [y0]])*exp(lambda*t)

>    evalm(A-lambda*I2)&*evalm(u0) = 0;

`&*`(matrix([[-1-lambda, 4], [-1, -1-lambda]]),matrix([[x0], [y0]])) = 0

>    chareqn :=
(A[1,1]-lambda)*(A[2,2]-lambda)-A[1,2]*A[2,1] = 0;
lam := solve(chareqn,lambda);

chareqn := (-1-lambda)^2+4 = 0

lam := -1+2*I, -1-2*I

Maple uses the symbol   I   for the complex number written elsewhere in Mathematics as i .
 The eigenvalues are both complex.


To find an eigenvector [x0,y0] corresponding to the first eigenvalue,
we consider the eigenproblem with lambda = lam[1]   :

>    lambda[1]:='lambda[1]': x0:='x0': y0:='y0':
evalm(A-lambda[1]*I2)&*evalm(u0) = 0;
evalm(A- lam[1]    *I2)&*evalm(u0) = 0;
evalm((A- lam[1]   *I2)&*u0) = 0;
eqn1:=(A[1,1]-lam[1])*x0 + A[1,2]*y0 = 0;
eqn2:= A[2,1]*x0 + (A[2,2]-lam[1])*y0 = 0;

`&*`(matrix([[-1-lambda[1], 4], [-1, -1-lambda[1]]]),matrix([[x0], [y0]])) = 0

`&*`(matrix([[-2*I, 4], [-1, -2*I]]),matrix([[x0], [y0]])) = 0

matrix([[-2*I*x0+4*y0], [-x0-2*I*y0]]) = 0

eqn1 := -2*I*x0+4*y0 = 0

eqn2 := -x0-2*I*y0 = 0

>    x0 := 1; y0:=solve(eqn1,y0); eval(eqn2);
v1 := colvector([x0,y0]); x0:='x0': y0:='y0':

x0 := 1

y0 := 1/2*I

0 = 0

v1 := matrix([[1], [1/2*I]])

>    t:='t':
u1 := v1*exp(lambda[1]*t);
u1 := evalm(v1)*exp(lam[1]*t);

u1 := v1*exp(lambda[1]*t)

u1 := matrix([[1], [1/2*I]])*exp((-1+2*I)*t)

>    e2It := cos(2*t)+I*sin(2*t):
u1 := evalm(v1)*exp(-t)*e2It;

u1 := matrix([[1], [1/2*I]])*exp(-t)*(cos(2*t)+sin(2*t)*I)

>    u1 := evalm(v1*exp(-t)*e2It);
u1:=map(f->combine(f),u1);

u1 := matrix([[exp(-t)*(cos(2*t)+sin(2*t)*I)], [1/2*I*exp(-t)*(cos(2*t)+sin(2*t)*I)]])

u1 := matrix([[exp(-t)*cos(2*t)+exp(-t)*sin(2*t)*I], [1/2*I*exp(-t)*cos(2*t)-1/2*exp(-t)*sin(2*t)]])

>    w1 := colvector([exp(-t)*cos(2*t),-(1/2)*exp(-t)*sin(2*t)]);
w2 := colvector([exp(-t)*sin(2*t),(1/2)*exp(-t)*cos(2*t)]);

w1 := matrix([[exp(-t)*cos(2*t)], [-1/2*exp(-t)*sin(2*t)]])

w2 := matrix([[exp(-t)*sin(2*t)], [1/2*exp(-t)*cos(2*t)]])

>    map(f->diff(f,t),w1) = evalm(A&*w1);
map(f->diff(f,t),w2) = evalm(A&*w2);

matrix([[-exp(-t)*cos(2*t)-2*exp(-t)*sin(2*t)], [1/2*exp(-t)*sin(2*t)-exp(-t)*cos(2*t)]]) = matrix([[-exp(-t)*cos(2*t)-2*exp(-t)*sin(2*t)], [1/2*exp(-t)*sin(2*t)-exp(-t)*cos(2*t)]])

matrix([[-exp(-t)*sin(2*t)+2*exp(-t)*cos(2*t)], [-1/2*exp(-t)*cos(2*t)-exp(-t)*sin(2*t)]]) = matrix([[-exp(-t)*sin(2*t)+2*exp(-t)*cos(2*t)], [-1/2*exp(-t)*cos(2*t)-exp(-t)*sin(2*t)]])

>    u_gen = c1*w1 + c2*w2;
u_gen = evalm(c1*w1 + c2*w2);

u_gen = c1*w1+c2*w2

u_gen = matrix([[c1*exp(-t)*cos(2*t)+c2*exp(-t)*sin(2*t)], [-1/2*c1*exp(-t)*sin(2*t)+1/2*c2*exp(-t)*cos(2*t)]])

where c1 and c2 are arbitrary real constants.

>    w1_x := w1[1,1]:  w2_x := w2[1,1]:
w1_y := w1[2,1]:  w2_y := w2[2,1]:
plot([w1_x,w2_x],t=0..10,labels=[`t`,`x`],title=`x-components of w1(red), w2(green)`);
plot([w1_y,w2_y],t=0..10,labels=[`t`,`y`],title=`y-components of w1(red) and of w2(green)`);

[Maple Plot]

[Maple Plot]

>    w1xyplot := plot([w1_x,w1_y,t=0..10],-1.5..1.5,-1..1,color=red):
w2xyplot := plot([w2_x,w2_y,t=0..10],color=green):
with(plots): display({w1xyplot,w2xyplot},labels=[`x`,`y`],title=`phase plane plots: w1(red), w2(green)`);

[Maple Plot]

>