option nolet
print"******************************************************"
print"* vector field 2d : with trajectory , version 1.0    *"
print"* peter.vlasschaert@gmail.com,10/10/2020             *"
print"******************************************************"
xa=5  ! horizontal,  values
ya=5  ! vertical  ,  values
! windows coordinate system : canvas
set window -xa,xa,-ya,ya
!parameters: ' draw vector field
w =0.05 ! vector field  - - - -, trajectory 'init to final'
a=2
b=-1
c=1
d=4
! draw axes sysem 'x-axes,y-axes'
plot -xa,0;xa,0 
plot 0,-ya;0,ya
!'=> trajectory parameters '
dim xt(0:1),yt(0:1)
nu_step = 1000   ! number of steps init value to end 
mat redim xt(0:nu_step+1),yt(0:nu_step+1)
!=> init value trajectory
xt(0) = 0.01 ! use 0 ,division by zero or -xa < xt(0)< xa
yt(0) = 0.01 ! use 0 ,division by zero or -ya < yt(0)< ya
set color 2 ! green
! '=> plot vector field 2D'
st = .25        ! real value
for x=-xa to xa step st
  for y=-ya to ya step st
  dx = a*x+b*y  ! ' vector field'
  dy = c*x+d*y
   ds=sqr(dx^2+dy^2)
 if ds >0 then plot x-w*dx/ds,y-w*dy/ds;x+w*dx/ds,y+w*dy/ds
  next y
next x
! '=> plot trajectory init -> final'
set color 4 ! red
for i=0 to nu_step !  'integer values'
    dx0 = a*xt(i)+b*yt(i)
    dy0 = c*xt(i)+d*yt(i)
    ds0 =sqr((dx0)^2+(dy0)^2)
       xt(i+1) = xt(i) + w*dx0/ds0
       yt(i+1) = yt(i) + w*dy0/ds0
    plot lines:xt(i),yt(i);xt(i+1),yt(i+1)
         xt(i) =xt(i+1)  ! connect lines,x values
         yt(i) =yt(i+1)  ! connect lines,y values
next i
end