---- Solving systems of two by two equations
-- By Mike Gage  gage@cc.rochester.edu  Apr 1995

-- The example solves an equation of the form y'=f(x,y) but the routines will solve any first order system of ODE's and plot the first two coordinates.  It's slow, but instructive since you can see all the workings and even change some of them to see what happens.

-- We define F given f(t,x)

F(T,Z)={1,f(Z[1],Z[2])}  -- non-parametric version

ic(x,y)=(stepto(x,y,X),{Z[1],Z[2]})


stepto(x,y,stop) = init(x,y) when stop=0,
            I:=0,(RungeKutta,  
            T:=T+dT,I:=I+1) while I<5

init(x,y)=Z:={x,y},T:=x,dT:=dt

Euler=Z:=Z+dT*F(T,Z)

RungeKutta=K1:=dT*F(T,Z),
 K2:=dT*F(T+dT/2,Z+K1/2),
 K3:=dT*F(T+dT/2,Z+K2/2),
 K4:=dT*F(T+dT,Z+K3),
 Z:=Z+(K1+2*K2+2*K3+K4)/6
-------------

axis={{Xmin,0},{Xmax,0},   -- t axis horizontal  
      {0/0,0/0},           -- lift pen
      {0,Ymin},{0,Ymax}}   -- x axis vertical

-- set window size
Xmin=-5; Xmax=5; Ymin=0; Ymax=2

-- set time step
dt=.05

-- set number of points to plot
Xsteps=20

--define direction field

f(t,x)=x-x^2

plotline axis

--plot integral curves  (ic)
plotline ic(-4, 0.01)
plotline ic(-3, 0.01)
plotline ic(-2, 0.01)
plotline ic(-1, 0.01)
plotline ic(0, 0.01)
plotline ic(1, 0.01)
plotline ic(-4,2)
plotline ic(-3,2)
plotline ic(-2,2)
plotline ic(-1,2)
plotline ic(0,2)