9.4.   Example 1:   x' = x - x^2 - xy, y' = 2y - y^2 - xy.  

>	with(DEtools):with(plots):with(plottools): # Load various packages

 

>	# First, I mark the three critical points with small blue disks.

 

>	C1:=disk([0,0],0.05,color=blue):

 

>	C2:=disk([0,2],0.05,color=blue):

 

>	C3:=disk([1,0],0.05,color=blue):

 

>	# Then I create the solutions for several initial conditions.

 

>

solutioncurves:=phaseportrait([D(x)(t) = x(t)-x(t)^2-x(t)*y(t), D(y)(t)= 2*y(t)-y(t)^2-x(t)*y(t)],[x(t),y(t)],t=0..3, [[x(0)=3,y(0)=3],[x(0)=2,y(0)=3],[x(0)=1,y(0)=3],[x(0)=0.5,y(0)=3],[x(0)=1,y(0)=0.1],[x(0)=0.5,y(0)=0.1],[x(0)=2,y(0)=0.1],[x(0)=0.1,y(0)=0.1],[x(0)=0.0,y(0)=0.1],[x(0)=3,y(0)=0],[x(0)=3,y(0)=1],[x(0)=3,y(0)=2]],x=0..3,y=0..3,linecolor=red,arrows=none,numpoints=100):

 

>	# Then I create the vector field.

 

>	vectorfield:=fieldplot([x-x^2-x*y,2*y-y^2-x*y],x=0..3,y=0..3,arrows=slim,anchor=tail,fieldstrength=maximal(3),grid=[20,20]):

 

>	# Finally, this is all displayed.

 

>	display(solutioncurves,vectorfield,C1,C2,C3);