#
# The user routines in this package are
#
#  fieldplot:	create a direction field plot for a first order ODE
#  phaseplot:	create a phase plot with flow lines for two first order ODE's
#  firsteuler:  solves a first order ODE using Eulers method
#  impeuler:	solves a first order ODE using the "improved Euler method"
#  rungekutta:	numerical solution to a system of first order ODEs using
#  		a 4'th order Runge-Kutta method.  The rungekuttahf routine
#		is a hardware floating point version of rungekutta
#
# There is a Maple V help file included for each of these routines.
# After reading in the procedures, to get help for phaseplot type ?phaseplot
# 
# June 1992.
# Fixed arrow routines which would not scale arrows in extreme cases.
# Changed (improved?) the shape of the arrows.
# Added type checking to improve the error messages the user would
# see when the parameters were incorrectly input, thanks to Michael 
# Monagan for his help on this.
# 
# November 1991.
# Added help files for rungekutta, rungekuttahf, firsteuler, impeuler.
# Fixed a bug in rungekuttahf which was restricting accuracy to current
# Digits setting, rather than the accuracy of evalhf.
#
# April 1991.
# This file contains Maple procedures which work with Maple V. They are
# used to analyze systems of first order differential equations.
# 
# Hard copy documentation for the procedures is available as chapter 5
# of a book published by Brooks/Cole Publishing Company called Maple V
# Flight Manual.
# 
# Thanks to George Labahn of the Waterloo group for sharing his code,
# parts of which helped to speed the following code. Also, the technique
# of putting arrows on vectors is his. Michael Monagan has been very
# helpful testing the routines and making suggestions for general
# improvements. Others in the Waterloo group have also been helpful for
# testing and suggestions.
# 
# Daniel Schwalbe
# 6635 County Road 101
# Corcoran, MN 55340
# 
# e-mail dschwalb@daisy.uwaterloo.ca
# 
#

`help/text/fieldplot` := TEXT( 
`FUNCTION: fieldplot - create a direction field with flow lines `,
`for first order differential equation.`,
``,
`CALLING SEQUENCE:`,
` fieldplot(f, h, v)`,
` fieldplot(f, h, v,{inits},...)`,
``,
`PARAMETERS:`,
` f - function of two variables where D(y)(t)=f(t,y).`,
` h - horizontal range`,
` v - vertical range`,
` {inits} - set of initial values for flow lines.`,
``,
`SYNOPSIS: - A typical call to the fieldplot function is`,
`fieldplot(f,a..b,c..d), where f is a real function of two`,
`variables and a..b,c..d specifies the horizontal and vertical real`,
`range. For each point, (tp,yp), of a grid a short line segment is`,
`drawn with slope f(tp,yp).`,
``,
`- Each initial point is specified as a list, [t0,y0]. A fourth`,
`order Runge/Kutta scheme is used to calculate the solution to`,
`D(y)(t)=f(t,y) through the point (t0,y0) over the entire`,
`horizontal range. In particular, the flow line will be`,
`erratic over any interval where a solution does not exist.In`,
`this case there is an option to use an adaptive numerical scheme`,
`which is much slower but will give a nice plot in some cases.`,
``,
`- Remaining arguments are interpreted as options which are`,
`specified as equations of the form option = value. There are`,
`currently 4 options -- grid, stepsize, iterations, and scheme.`,
`The value for grid is a list, [m,n], where the current default is`,
`[16,12]. The value for stepsize controls the step size in the`,
`Runge/Kutta numerical scheme and determines how many points will`,
`be plotted, the current default is (b-a)/20 and 20 points are`,
`plotted. It is possible to decrease the stepsize without`,
`increasing the number of points plotted (which slows down the`,
`routine and leads to unreasonably large plot structures) by`,
`setting the value for iterations. iterations = 10 decreases the`,
`stepsize by a factor of 10, and every 10th point will be stored and`,
`plotted. The option scheme=adaptive specifies that an adaptive`,
`fourth-fifth Runge/Kutta numerical scheme will be used. This `,
`should be used carefully because it will be much slower.`,
``,
``,
`EXAMPLES: fieldplot((t,y)->t-y,-4..4,-3..3);`,
`fieldplot((t,y)->y,0..1,0..3,{[0,1]});`,
`fieldplot((t,y)->-t/y,-4..4,-4..4,{[0,1],[0,2],[0,3],[0,4]});`,
`fieldplot((t,y)->-t/y,-4..4,-4..4,[0,1],scheme=adaptive);`,
`fieldplot((t,y)->sin(t+y),-4..4,-4..4,grid=[16,14]);`,
``,
`eq:=(t,y)->sin(t+y):`,
`inits:={seq(seq([2*i,2*j],i=-2..2),j=-2..2)}:`,
`fieldplot(eq,-4..4,-4..4,inits,grid=[0,0]);`,
``,
`fieldplot((t,y)->-y^2+t,0..20,0..10,[0,1]);`,
`fieldplot((t,y)->-y^2+t,0..20,0..10,[0,1],iterations=10);`,
`fieldplot((t,y)->-y^2+t,0..20,0..10,[0,1],scheme=adaptive);`,
``,
`eq4:=(x,y)->(1-x^2-y^2)/(y-x+2):`,
`fieldplot(eq4,-2..2,-2..2,grid=[15,15]);`,
``,
`eq4:=(x,y)->(1-x^2-y^2)/(y-x+2):`,
`init4:=seq(seq([j,i],i=-2..2),j=-2..2):`,
`fieldplot(eq4,-2..2,-2..2,{init4});`,
``,
`eq4:=(x,y)->(1-x^2-y^2)/(y-x+2):`,
`fieldplot(eq4,-2..2,-2..2,[1,1],scheme=adaptive);`,
``,
`eq5:=(t,y)->y*tan(t);`,
`init5:=seq(seq([i,j],i=-2..2),j=-2..2):`,
`fieldplot(eq5,-6..6,-6..6,grid=[0,0],{init5},`,
` iterations=5,stepsize=0.1);`,
``,
`eq6:=(t,y)->-y*tan(t);`,
`init6:=seq(seq([2*i,2*j],i=-2..2),j=-2..2):`,
`fieldplot(eq6,-6..6,-6..6,grid=[0,0],{init6},`,
` iterations=5,stepsize=.1);`,
` `,
`see also: firsteuler, impeuler, rungekutta, rungekuttahf, phaseplot`):

`help/text/phaseplot` := TEXT( 
`FUNCTION: phaseplot - draw phase space with flow lines for second`,
` order autonomous differential equation.`,
` `,
`CALLING SEQUENCE:`,
` phaseplot([f1,f2], h, v)`,
` phaseplot([f1,f2], h, v,{inits},...)`,
``,
`PARAMETERS: f1,f2 - two functions of three variables where`,
`D(x)(t)=f1(t,x,y) and D(y)(t)=f2(t,x,y) where normally f1 and f2`,
`do not depend on t. t=0 is substituted to draw the phase space`,
`otherwise. Note it will still be interesting to draw flow`,
`lines in the non-autonomous case, i.e. where f1 and f2 depend`,
`explicity on time. h - horizontal range v - vertical range {inits}`,
`- set of initial values for flow lines.`,
``,
`SYNOPSIS: - A typical call to the phaseplot function is`,
`phaseplot([f1,f2],a..b,c..d), where f1 and f2 are real functions of`,
`three variables and a..b,c..d specifies the horizontal and`,
`vertical real range. For each point, (xp,yp), of a grid, a scaled`,
`vector with basepoint (xp,yp) and endpoint`,
`(f1(0,xp,yp),f2(0,xp,yp)) is drawn.`,
``,
`- Each initial point is specified as a list, [t0,x0,y0]. A`,
`fourth order Runge/Kutta scheme is used to calculate the solution`,
`to the system of equations, D(x)(t)=f1(t,x,y),`,
`D(y)(t)=f2(t,x,y),x(t0)=x0,y(t0)=y0. The default is to use a step`,
`size of 0.1 and do 50 iterations.`,
``,
`- Remaining arguments are interpreted as options which are`,
`specified as equations of the form option = value. There are`,
`currently 4 options -- grid, stepsize, numsteps, and iterations.`,
`The value for grid is a list, [m,n], where the current default is`,
`[16,14]. The value for stepsize controls the step size in the`,
`Runge/Kutta numerical scheme, the current default is 0.1. The`,
`value for numsteps specifies how many steps forward and backward`,
` will be calculated and stored by the Runge/Kutta routine. The`,
`current default is 50 forward and back. It is possible to`,
`decrease the stepsize without increasing the number of points`,
`plotted (which slows down the routine and leads to unreasonably`,
`large plot structures) by setting the value for iterations.`,
`iterations = 10, would decrease the stepsize by a factor of 10`,
`and every 10th point would be stored and plotted.`,
` `,
``,
`EXAMPLES: `,
`phaseplot([(t,x,y)->-x,(t,x,y)->y],-4..4,-4..4);`,
``,
`eq1:=(t,x,y)-> 5*x-x^2-x*y: `,
`eq2:=(t,x,y)->x*y-2*y:`,
`inits:={seq([0,5,i],i=1..9),seq([0,i,0.1],i=1..9)}:`,
`phaseplot([eq1,eq2],0..10,0..10,inits,grid=[15,15]);`,
``,
`readlib(Heaviside): `,
`eq3:=(t,y,v)->v: `,
`eq4:=(t,y,v)->-Heaviside(y):`,
`phaseplot([eq3,eq4],-4..4,-4..4); `,
``,
`readlib(Heaviside): `,
`eq5:=(t,y,v)->v:`,
`eq6:=(t,y,v)->-Heaviside(y-5)*v/2-(1-Heaviside(y-5))*v/10+1:`,
`phaseplot([eq5,eq6],0..10,0..5,[0,0,0],iterations=10,`,
` numsteps=100,grid=[15,15]);`,
``,
`eq7:=(t,x,y)-> -x+x*y: `,
`eq8:=(t,x,y)->-y+2*x*y:`,
`inits2:={seq(seq([0,i/2,j/2],i=0..2),j=0..2)}:`,
`phaseplot([eq7,eq8],-.5..1.5,-.5..1.5,inits2,grid=[15,15]);`,
``,
`pend3:=(t,y,v) -> v;`,
`pend4:=(t,y,v) -> -sin(y) -0.2*v;`,
`pendinits:=seq([0,j/2],j=2..8):`,
`phaseplot([pend3,pend4],-10..10,-4..4,{pendinits},`,
` numsteps=200);`,
``,
`ropedv:= proc(t,v,y) if v > 0 then (9800 -3*980*y - 3*v^2)/(3*y)`,
`   else (9800 -3*980*y)/(3*y) fi end;`,
`ropedy:= (t,v,y) -> v;`,
`phaseplot([ropedv,ropedy],-59..59,0.1..5);`,
``,
`eq11:=(t,y,v)->v^2-y^2;`,
`eq12:=(t,y,v)->y-2*v;`,
`eq1init:=seq(seq([i,j],j=-2..2),i=-2..2):`,
`phaseplot([eq11,eq12],-5..5,-5..5,grid=[12,12],{eq1init});`,
``,
`eq21:=(t,x,y)->cos(x+y);`,
`eq22:=(t,x,y)->-sin(y);`,
`phaseplot([eq21,eq22],-6..6,-6..6,grid=[25,25]);`,
``,
`see also: firsteuler, impeuler, rungekutta, rungekuttahf, fieldplot`
):

`help/text/rungekutta` := TEXT( 
`FUNCTION: rungekutta - approximate a solution to a system of first`,
`order ordinary differential equations using a fourth order Runge/Kutta`,
`numerical method. Solution is returned as an array of lists.`,
`rungekuttahf is an alternative form of this procedure which uses`,
`hardware floating point operations and is not affected by the Digits`,
`settings. rungekuttahf is limited to double precision on most systems`,
`and will be about as accurate as rungekutta using a Digits setting `,
`of 15. It is however, MUCH faster.`,
``,
`CALLING SEQUENCE:`,
` rungekutta(fcnlist,initvalue,stepsize,numsteps)`,
` rungekutta(fcnlist,initvalue,stepsize,numsteps,iterations)`,
``,
`PARAMETERS:`,
` fcnlist - list of functions, [f1,f2,...,fn], of n+1 real variables.  `,
`           The differential equation to be solved is of the form:`,
`				dx1/dt = f1(t,x1,x2,...,xn)`,
`				dx2/dt = f2(t,x1,x2,...,xn)`,
`				   .   .          .`,
`				   .   .          .`,
`				   .   .          .`,
`				dxn/dt = fn(t,x1,x2,...,xn)`,
` initvalue - list of n+1 numbers, [i0,i1,i2,...,in],  which represents `,
`             the initial value of the differential equation. `,
` stepsize - step size used by the numerical method.`,
` numsteps - number of steps taken by numerical method.`,
` iterations - optional value which decreases the step size without `,
`              storing more points.  New step size will be `,
`			  stepsize/iterations.`,
``,
`SYNOPSIS: - A typical call to the rungekutta function is`,
`rungekutta([f1,f2,...,fn],[i0,i1,i2,...,in],stepsize,numsteps), where`,
`f1,f2,...,fn are real valued functions of n+1 variables, stepsize is`,
`the step size used by the numerical method, and numsteps specifies the`,
`number of steps to be taken by the numerical method. An array of (n+1)`,
`lists of (numsteps+1) numbers is returned, where each list represents`,
`a point calculated by the numerical method.`,
``,
`- An optional fifth argument is used to reduce the step size without `,
`storing more points.  This is especially useful when graphing `,
`solutions which require a small step size and a large plot structure `,
`is undesirable.  Also useful when only the final value is wanted.`,
``,
`- makelist is a procedure in the ODE file to extract a list of values`,
`from an array of lists. This is useful for making two dimensional`,
`plots of solutions calculated by rungekutta.  makelist(A,m,n) will `,
`make a list of values from the mth and nth components from each `,
`entry of A. m defaults to 1 and n defaults to 2 if no values are `,
`entered.`,
``,
`EXAMPLES: `,
`rkpts1:=rungekutta((t,y)->t-y,[0,1],.1,10);`,
`plot(makelist(rkpts1));`,
``,
`eq1:=(t,x,y)-> y;`,
`eq2:=(t,x,y)->-x;`,
`rkpts2:=rungekutta([eq1,eq2],[0,1,1],.1,60);`,
`plot(makelist(rkpts2,2,3));`,
``,
`dp1:=(t,x,xp,y,yp)->xp;`,
`dp2:=(t,x,xp,y,yp)->-2*sin(x)+sin(y);`,
`dp3:=(t,x,xp,y,yp)->yp;`,
`dp4:=(t,x,xp,y,yp)->-2*sin(y)+2*sin(x);`,
`dpinit:=[0,0,0,0,2];`,
`dppts:=rungekutta([dp1,dp2,dp3,dp4],dpinit,1,10,1):`,
`dppts2:=rungekutta([dp1,dp2,dp3,dp4],dpinit,1,10,2):`,
`dppts3:=rungekutta([dp1,dp2,dp3,dp4],dpinit,1,10,4):`,
`plot({makelist(dppts,2,4),makelist(dppts2,2,4),`,
`  makelist(dppts3,2,4)});`,
``,
`see also: firsteuler, impeuler, rungekuttahf, fieldplot, phaseplot`):

`help/text/rungekuttahf` := TEXT( 
`FUNCTION: rungekuttahf - approximate a solution to a system of up to`,
`four first order ordinary differential equations using a fourth order`,
`Runge/Kutta numerical method. Solution is returned as an array of`,
`lists. rungekutta is an alternative form of this procedure which uses`,
`Maple's arbritrary precision floating point operations and is thus`,
`affected by the Digits settings. rungekutta is useful for studying the`,
`effect of the round off error that occurs with numerical methods or`,
`with differential equations for which rungekuttahf fails. rungekuttahf`,
`is limited to double precision on most systems and will be about as`,
`accurate as rungekutta using a Digits setting of 15. rungekuttahf is`,
`however, MUCH faster.`,
``,
`CALLING SEQUENCE:`,
` rungekuttahf(fcnlist,initvalue,stepsize,numsteps)`,
` rungekuttahf(fcnlist,initvalue,stepsize,numsteps,iterations)`,
``,
`PARAMETERS:`,
` fcnlist - list of functions, [f1,f2,...,fn], of n+1 real variables.  `,
`           The differential equation to be solved is of the form:`,
`				dx1/dt = f1(t,x1,x2,...,xn)`,
`				dx2/dt = f2(t,x1,x2,...,xn)`,
`				   .   .          .`,
`				   .   .          .`,
`				   .   .          .`,
`				dxn/dt = fn(t,x1,x2,...,xn)`,
` initvalue - list of n+1 numbers, [i0,i1,i2,...,in],  which represents `,
`             the initial value of the differential equation. `,
` stepsize - step size used by the numerical method.`,
` numsteps - number of steps taken by numerical method.`,
` iterations - optional value which decreases the step size without `,
`              storing more points.  New step size will be `,
`			  stepsize/iterations.`,
``,
`SYNOPSIS: - A typical call to the rungekuttahf function is`,
`rungekuttahf([f1,f2,...,fn],[i0,i1,i2,...,in],stepsize,numsteps), where`,
`f1,f2,...,fn are real valued functions of n+1 variables, stepsize is`,
`the step size used by the numerical method, and numsteps specifies the`,
`number of steps to be taken by the numerical method. An array of (n+1)`,
`lists of (numsteps+1) numbers is returned, where each list represents`,
`a point calculated by the numerical method.`,
``,
`- An optional fifth argument is used to reduce the step size without `,
`storing more points.  This is especially useful when graphing `,
`solutions which require a small step size and a large plot structure `,
`is undesirable.  Also useful when only the final value is wanted.`,
``,
`- makelist is a procedure in the ODE file to extract a list of values`,
`from an array of lists. This is useful for making two dimensional`,
`plots of solutions calculated by rungekuttahf.  makelist(A,m,n) will `,
`make a list of values from the mth and nth components from each `,
`entry of A. m defaults to 1 and n defaults to 2 if no values are `,
`entered.`,
``,
`EXAMPLES: `,
`rkpts1:=rungekuttahf((t,y)->t-y,[0,1],.1,10);`,
`plot(makelist(rkpts1));`,
``,
`eq1:=(t,x,y)-> y;`,
`eq2:=(t,x,y)->-x;`,
`rkpts2:=rungekuttahf([eq1,eq2],[0,1,1],.1,60);`,
`plot(makelist(rkpts2,2,3));`,
``,
`dp1:=(t,x,xp,y,yp)->xp;`,
`dp2:=(t,x,xp,y,yp)->-2*sin(x)+sin(y);`,
`dp3:=(t,x,xp,y,yp)->yp;`,
`dp4:=(t,x,xp,y,yp)->-2*sin(y)+2*sin(x);`,
`dpinit:=[0,0,0,0,2];`,
`dppts:=rungekuttahf([dp1,dp2,dp3,dp4],dpinit,.1,100):`,
`plot({makelist(dppts,2,4)});`,
``,
`eq3:=(t,x,y)->x*sin(t)+5*y;`,
`eq4:=(t,x,y)->-10*x+y*sin(t);`,
`rkpts2:=rungekuttahf([eq3,eq4],[0,2,0],.1,100,5);`,
`plot(makelist(rkpts2,2,3));`,
`plot({makelist(rkpts2,1,2),makelist(rkpts2,1,3)});`,
`#for the Maple expert and the brave alike, construct a PLOT3D structure`,
`plot1:=PLOT3D(CURVELIST([seq([rkpts2[i][1],rkpts2[i][2],rkpts2[i][3]],`,
`    i=0..100)]),AXES(BOXED)):`,
`plot1;`,
`	`,
``,
`see also: firsteuler, impeuler, rungekutta, fieldplot, phaseplot`):

`help/text/firsteuler` := TEXT( 
`FUNCTION: firsteuler - approximate a solution to a first order`,
`ordinary differential equations using Euler's method.`,
`Solution is returned as an array of lists. `,
``,
`CALLING SEQUENCE:`,
` firsteuler(f,i,h,n)`,
``,
`PARAMETERS:`,
` f - function of 2 real variables.  The differential equation to be `,
`     solved is dy/dt = f(t,y)  `,
` i - list of 2 numbers. `,
` h - step size used by the numerical method.`,
` n - number of steps taken by numerical method.`,
``,
`SYNOPSIS: - A typical call to the firsteuler function is`,
`firsteuler(f,i,h,n) where f is a real valued functions of 2 variables,`,
`h is the step size used by the numerical method, and n specifies the`,
`number of steps to be taken by the numerical method. An array of (n+1)`,
`lists of 2 numbers is returned, where each list represents a point`,
`calculated by the numerical method.`,
``,
`- makelist is a procedure in the ODE file to extract a list of values`,
`from an array of lists. This is useful for making a plot of the`,
`approximate solution calculated by firsteuler. makelist(A) will make a`,
`list of values from the entries of A. `,
``,
`EXAMPLES: `,
`eulerpts:=firsteuler((t,y)->t-y,[0,1],.1,10);`,
`plot(makelist(eulerpts));`,
``,
`eq:=(t,y)->y;`,
`pts1:=firsteuler(eq,[0,1],.5,2);`,
`pts2:=firsteuler(eq,[0,1],.25,4);`,
`pts3:=firsteuler(eq,[0,1],.1,10);`,
`plot({makelist(pts1),makelist(pts2),makelist(pts3),`,
`      exp(t)},0..1,style=LINE);`,
`	  `,
`eq2:=(t,y)->-t/y;`,
`eq2pts:=firsteuler(eq2,[0,4],.2,25):`,
`plot(makelist(eq2pts));`,
``,
`eq3:=(t,y)->5*y-6*exp(-t);`,
`eq3pts:=firsteuler(eq3,[0,1],.1,10);`,
`eq3sol:=dsolve({diff(y(t),t)=eq3(t,y),y(0)=1},y(t));`,
`plot({makelist(eq3pts),rhs(eq3sol)},0..1);`,
``,
`eq4:=(t,y)->25*y*(1-y);`,
`eq4pts:=firsteuler(eq4,[0,1.3],.1,10):`,
`plot(makelist(eq4pts));`,
``,
`eq4:=(t,y)->25*y*(1-y);`,
`eq4pts:=firsteuler(eq4,[0,1.3],.1,10):`,
`eq4pts2:=firsteuler(eq4,[0,1.3],.05,20):`,
`eq4pts3:=firsteuler(eq4,[0,1.3],.025,40):`,
`plot({makelist(eq4pts),makelist(eq4pts2),makelist(eq4pts3)});`,
``,
``,
`see also: impeuler, rungekutta, rungekuttahf, fieldplot, phaseplot`):

`help/text/impeuler` := TEXT(
`FUNCTION: impeuler - approximate a solution to a first order ordinary`,
`differential equations using the improved Euler method. Solution is`,
`returned as an array of lists. `,
``,
`CALLING SEQUENCE:`,
` impeuler(f,i,h,n)`,
``,
`PARAMETERS:`,
` f - function of 2 real variables.  The differential equation to be `,
`     solved is dy/dt = f(t,y)  `,
` i - list of 2 numbers. `,
` h - step size used by the numerical method.`,
` n - number of steps taken by numerical method.`,
``,
`SYNOPSIS: - A typical call to the impeuler function is`,
`impeuler(f,i,h,n) where f is a real valued functions of 2 variables,`,
`h is the step size used by the numerical method, and n specifies the`,
`number of steps to be taken by the numerical method. An array of (n+1)`,
`lists of 2 numbers is returned, where each list represents a point`,
`calculated by the numerical method.`,
``,
`- makelist is a procedure in the ODE file to extract a list of values`,
`from an array of lists. This is useful for making a plot of the`,
`approximate solution calculated by impeuler. makelist(A) will make a`,
`list of values from the entries of A. `,
``,
`EXAMPLES: `,
`impts:=impeuler((t,y)->t-y,[0,1],.1,10);`,
`plot(makelist(impts));`,
``,
`eq:=(t,y)->y;`,
`impts1:=impeuler(eq,[0,1],.5,2);`,
`impts2:=impeuler(eq,[0,1],.25,4);`,
`impts3:=impeuler(eq,[0,1],.1,10);`,
`plot({makelist(impts1),makelist(impts2),makelist(impts3),`,
`      exp(t)},0..1,style=LINE);`,
`	  `,
`eq2:=(t,y)->-t/y;`,
`eq2impts:=impeuler(eq2,[0,4],.2,25):`,
`plot(makelist(eq2impts));`,
``,
`eq3:=(t,y)->5*y-6*exp(-t);`,
`eq3impts:=impeuler(eq3,[0,1],.1,10);`,
`eq3sol:=dsolve({diff(y(t),t)=eq3(t,y),y(0)=1},y(t));`,
`plot({makelist(eq3impts),rhs(eq3sol)},0..1);`,
``,
`eq4:=(t,y)->25*y*(1-y);`,
`eq4impts:=impeuler(eq4,[0,1.3],.1,10):`,
`plot(makelist(eq4impts));`,
``,
`eq4:=(t,y)->25*y*(1-y);`,
`eq4impts:=impeuler(eq4,[0,1.3],.1,10):`,
`eq4impts2:=impeuler(eq4,[0,1.3],.05,20):`,
`eq4impts3:=impeuler(eq4,[0,1.3],.025,40):`,
`eq4implot:=plot({makelist(eq4impts),makelist(eq4impts2),`,
`              makelist(eq4impts3)}):`,
`eq4field:=fieldplot(eq4,0..1,0.8..1.5):`,
`with(plots):`,
`display({eq4implot,eq4field});`,
``,
``,
`see also: firsteuler, rungekutta, rungekuttahf, fieldplot, phaseplot`):


#################################################

firsteuler := proc(f,initlist,h0,n,iterations)
local i,xk,yk,pts,h,npts;
options `Copyright 1991 by Daniel Schwalbe`;
pts:=array(0..n);
if nargs>4 then
   npts:=iterations;
else npts:=1 fi;
h:=evalf(h0/npts);  
xk:=evalf(initlist[1]); yk:=evalf(initlist[2]);
pts[0] := [xk,yk];
for i to n while type(yk,numeric) do
  to npts while type(yk,numeric) do
    yk:=yk + h*f(xk,yk);
    xk:=xk+h od;
  pts[i]:=[xk,yk] od:
if i <= n then ERROR(`undefined procedure`)
else op(pts) fi:
end:

##################################################

impeuler:=proc(f,initlist,h0,n,iterations)
local h,xk,yk,pts,j,npts:
options `Copyright 1991 by Daniel Schwalbe`;
pts:=array(0..n);
if nargs>4 then 
  npts:=iterations
  else npts:=1 fi;
h:=evalf(h0/npts):
xk:=evalf(initlist[1]); yk:=evalf(initlist[2]);
pts[0] := [xk,yk];
for j to n while type(yk,numeric) do
  to npts while type(yk,numeric) do
    yk:=yk+h/2*(f(xk,yk)+f(xk+h,yk+h*f(xk,yk))):
    xk:=xk+h od:
  pts[j]:=[xk,yk] od:
if j <= n then ERROR(`undefined procedure`)
 else op(pts) fi:
end:
##################################################
rungekutta:=proc(f,initlist,h0,n,iterations)
local k1,k2,k3,k4,fl,h,pts,j,i,k,m,npts,temp:
options `Copyright 1991 by Daniel Schwalbe`;
m:=nops([op(initlist)])-1;
fl:=f; fl:=[op(fl)]; pts:=array(0..n); 
if nargs > 4 then
  npts:=iterations
  else npts:=1 fi;
h:=evalf(h0/npts): 
pts[0] := map(evalf,[op(initlist)]);
temp:=pts[0];
for j from 0 to n-1 while 
  not member(false,map(type,(temp,numeric))) do
    to npts while not member(false,map(type,(temp,numeric))) do
      for i to m do k1[i]:=
         evalf(fl[i](op(temp))) od:
      for i to m do k2[i]:= evalf(fl[i](temp[1]+h/2,
         seq(temp[k+1] + h*k1[k]/2,k=1..m))) od:
      for i to m do k3[i]:=evalf(fl[i](temp[1]+h/2,
         seq(temp[k+1] + h*k2[k]/2,k=1..m))) od:
      for i to m do k4[i]:=evalf(fl[i](pts[j][1]+h,
         seq(temp[k+1] + h*k3[k],k=1..m))) od:
    temp:=[temp[1]+h, seq(temp[i+1]+
      h*(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6,i=1..m)] od:
  pts[j+1]:= temp od:
if j<=n-1 then ERROR(`undefined procedure`) 
 else op(pts) fi:
end:

##################################################
rkhf:=proc(Akf,B,h,n,npts)
local j,i,AA,k:
options `Copyright 1991 by Daniel Schwalbe`;
AA:=array(0..npts,1..2);
for i to 2 do AA[0,i]:=B[0,i] od;
for i to n do
  for j to npts  do
     AA[j,1]:=AA[j-1,1]+h;
     AA[j,2]:=Akf(AA,j-1) od;
  for k to 2 do B[i,k]:=AA[npts,k] od;
  for k to 2 do AA[0,k]:=AA[npts,k] od;
od;
end:

rkhf2:=proc(Akf,Akf2,B,h,n,npts)
local j,i,AA,k:
options `Copyright 1991 by Daniel Schwalbe`;
AA:=array(0..npts,1..3);
for i to 3 do AA[0,i]:=B[0,i] od;
for i to n do
  for j to npts  do
     AA[j,1]:=AA[j-1,1]+h;
     AA[j,2]:=Akf(AA,j-1);
     AA[j,3]:=Akf2(AA,j-1) od:
  for k to 3 do B[i,k]:=AA[npts,k] od;
  for k to 3 do AA[0,k]:=AA[npts,k] od;
od;
end:

rkhf3:=proc(Akf,Akf2,Akf3,B,h,n,npts)
local j,i,AA,k:
options `Copyright 1991 by Daniel Schwalbe`;
AA:=array(0..npts,1..4);
for i to 4 do AA[0,i]:=B[0,i] od;
for i to n do
  for j to npts  do
     AA[j,1]:=AA[j-1,1]+h;
     AA[j,2]:=Akf(AA,j-1);
     AA[j,3]:=Akf2(AA,j-1):
     AA[j,4]:=Akf3(AA,j-1) od:
  for k to 4 do B[i,k]:=AA[npts,k] od;
  for k to 4 do AA[0,k]:=AA[npts,k] od;
od;
end:

rkhf4:=proc(Akf,Akf2,Akf3,Akf4,B,h,n,npts) 
local j,i,AA,k:
options `Copyright 1991 by Daniel Schwalbe`;
AA:=array(0..npts,1..5);
for i to 5 do AA[0,i]:=B[0,i] od;
for i to n do
  for j to npts  do
     AA[j,1]:=AA[j-1,1]+h;
     AA[j,2]:=Akf(AA,j-1);
     AA[j,3]:=Akf2(AA,j-1):
     AA[j,4]:=Akf3(AA,j-1):
     AA[j,5]:=Akf4(AA,j-1) od:
  for k to 5 do B[i,k]:=AA[npts,k] od;
  for k to 5 do AA[0,k]:=AA[npts,k] od;
od;
end:

rungekuttahf:=proc(f,init,h0,n,iterations)
local h,i,m,A,B,fl,k,k1,k2,k3,k4,precision,
    npts,Ak,Akf,Akf2,Akf3,Akf4,j;
options `Copyright 1991 by Daniel Schwalbe`;
precision := max(Digits,trunc(evalhf(Digits)));
m:=nops([op(init)])-1;
A:=array(0..n,1..m+1); 
if nargs > 4 then npts:=iterations
  else npts:=1 fi;
h:=evalf(h0/npts,precision):
for i to m+1 do A[0,i]:=evalf(init[i],precision) od;
if type(f,list) then fl:=f else fl:=[f] fi;
for i to m do k1[i]:=fl[i](seq(A[j,k],k=1..m+1)) od:
for i to m do k2[i]:=fl[i](A[j,1]+h/2,
   seq(A[j,k+1] + h*k1[k]/2,k=1..m)) od:
for i to m do k3[i]:=fl[i](A[j,1]+h/2,
   seq(A[j,k+1] + h*k2[k]/2,k=1..m)) od:
for i to m do k4[i]:=fl[i](A[j,1]+h,
   seq(A[j,k+1] + h*k3[k],k=1..m)) od:
for i to m do Ak[i]:=
   A[j,i+1]+h*(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6 od:
if m=1 then
   Akf:=unapply(Ak[1],A,j);
   evalhf(rkhf(Akf,var(A),h,n,npts))
elif m=2 then
   Akf:=unapply(Ak[1],A,j);
   Akf2:=unapply(Ak[2],A,j);
   evalhf(rkhf2(Akf,Akf2,var(A),h,n,npts))
elif m=3 then
   Akf:=unapply(Ak[1],A,j);
   Akf2:=unapply(Ak[2],A,j);
   Akf3:=unapply(Ak[3],A,j);
   evalhf(rkhf3(Akf,Akf2,Akf3,var(A),h,n,npts))
elif m=4 then
   Akf:=unapply(Ak[1],A,j);
   Akf2:=unapply(Ak[2],A,j);
   Akf3:=unapply(Ak[3],A,j);
   Akf4:=unapply(Ak[4],A,j);
   evalhf(rkhf4(Akf,Akf2,Akf3,Akf4,var(A),h,n,npts))
else print(`rungekuttahf is not designed for more than 4 equations.`);
     print(`Try rungekutta as an alternate.`) fi:
B:=array(0..n);
for i from 0 to n do B[i]:=[seq(A[i,k],k=1..m+1)] od;
op(B)
end:

##################################################

makelist:=proc(A,c10,c20)
local c1,c2,ri,rf,i;
options `Copyright 1991 by Daniel Schwalbe`;
if nargs = 1 then
   ri:=op(1,op(2,eval(A)));
   rf:=op(2,op(2,eval(A)));
   c1:=1;
   c2:=2;
else
   c1:=c10;
   c2:=c20;
   ri:=op(1,op(2,eval(A)));
   rf:=op(2,op(2,eval(A))) fi;
map(op,[seq([A[i][c1],A[i][c2]],i=ri..rf)]);
end:

#################################################
intbox:=proc(p,q,r,s,pta,pt)
local xnew,ynew;
options `Copyright 1991 by Daniel Schwalbe`;
if pt[1]<p then 
  xnew:=p;
  ynew:=(pt[2]-pta[2])/(pt[1]-pta[1])*(xnew-pta[1])+pta[2];
  pt[1]:=xnew; pt[2]:=ynew 
elif pt[1]>q then
  xnew:=q;
  ynew:=(pt[2]-pta[2])/(pt[1]-pta[1])*(xnew-pta[1])+pta[2];
  pt[1]:=xnew; pt[2]:=ynew fi;
if pt[2]<r then
  ynew:=r;
  xnew:=(ynew-pta[2])*(pt[1]-pta[1])/(pt[2]-pta[2])+pta[1];
  pt[1]:=xnew; pt[2]:=ynew 
elif pt[2]>s then
  ynew:=s;
  xnew:=(ynew-pta[2])*(pt[1]-pta[1])/(pt[2]-pta[2])+pta[1];
  pt[1]:=xnew; pt[2]:=ynew fi;
end:

##################################################
rkstep:=proc(f,hp,endpt,npts,pt)
local h,t,y,k1,k2,k3,k4;
options `Copyright 1991 by Daniel Schwalbe`;
t:=pt[1]:y:=pt[2]: h:=hp:
to npts while ((h>0 and t<endpt) or (h<0 and t>endpt)) do
  if (h>0 and t+h>endpt) or (h<0 and t+h<endpt) then
    h:=endpt-t fi:
    k1:=f(t,y):
    k2:=f(t+h/2,y+ h*k1/2):
    k3:=f(t+h/2,y+ h*k2/2):
    k4:=f(t+h,y+ h*k3):
    y:=y+h*(k1+2*k2+2*k3+k4)/6:
    t:=t+h: od:
pt[1]:=t: pt[2]:=y:
end:

fieldrk:=proc(A,f,CC,npts)
local pta,pt,j,h,p,q,r,s;
options `Copyright 1991 by Daniel Schwalbe`;
pt:=array(1..2):
pta:=array(1..2):
h:=CC[3]/npts;
r:=CC[6];s:=CC[7];
p:=CC[4];q:=CC[5];
pt[1]:=CC[1]: pt[2]:=CC[2]:
A[0]:=CC[1]; A[1]:=CC[2];
j:=0;CC[9]:=1:
while pt[1] < q and pt[2]<=s and pt[2]>=r do
  A[2*j+1]:=pt[2]: A[2*j]:=pt[1];
  CC[9]:=2*j+1:
  j:=j+1;
  rkstep(f,h,q,npts,pt): od:
if j>0 then
  pta[1]:=A[2*j-2];
  pta[2]:=A[2*j-1];
  intbox(p,q,r,s,pta,pt);
  A[2*j]:=pt[1];
  A[2*j+1]:=pt[2];
  CC[9]:=2*j+1 fi;
j:=0;CC[8]:=0;
pt[1]:=CC[1]; pt[2]:=CC[2];
while pt[1] > p and pt[2]>=r and pt[2]<=s do
  A[2*j+1]:=pt[2]: A[2*j]:=pt[1];
  CC[8]:=2*j:
  j:=j-1;
  rkstep(f,-h,p,npts,pt): od:
  if j<0 then
    pta[1]:=A[2*j+2];
    pta[2]:=A[2*j+3];
    intbox(p,q,r,s,pta,pt);
    A[2*j]:=pt[1];
    A[2*j+1]:=pt[2]; 
    CC[8]:=2*j fi;
end:

adaptiverk:=proc(f,CC)
local F,j,tt,Fval,yy,A,p,q,h,t,y;
options `Copyright 1991 by Daniel Schwalbe`;
p:=CC[4]: q:=CC[5]:
tt:=CC[1]:yy:=CC[2]:h:=CC[3]:
A[0]:=tt,yy:
F:=dsolve({D(y)(t)=f(t,y(t)),y(tt)=yy},y(t),numeric);
j:=0; Fval:=traperror(F(tt+h));
while Fval<>lasterror and tt+h<=q do
  tt:=tt+h; j:=j+1;
  A[j]:=Fval:
  Fval:=traperror(F(tt+h)) od;
CC[9]:=j:
j:=0;tt:=CC[1]: Fval:=traperror(F(tt-h));
while Fval<>lasterror and tt-h>=p do
  tt:=tt-h; j:=j-1;
  A[j]:=Fval:
  Fval:=traperror(F(tt-h)) od;
CC[8]:=j:
op(A)
end:

computelineshf:=proc(f,Ax,Ay,CC)
local m,n,dx,dy,k,slp,xinc,yinc,a,b,p,q,r,s,i,j,ldx,ldy;
m:=CC[1]; n:=CC[2];
p:=CC[4]; q:=CC[5];
r:=CC[6]; s:=CC[7];
dx:=(q-p)/m; dy:=(s-r)/n;
ldx:=3*dx/8; ldy:=3*dy/8;
a:=p;
for i from 0 to m do
  b:=r;
  for j from 0 to n do
    k:=2*j;
    slp:=f(a,b);
    if slp > dy/dx or -slp > dy/dx then
      xinc:=ldy/slp; yinc:=ldy 
    else xinc:=ldx; yinc:=ldx*slp fi;
    Ax[i,k]:=a-xinc;
    Ay[i,k]:=b-yinc;
    Ax[i,k+1]:=a+xinc;
    Ay[i,k+1]:=b+yinc;
    b:=b+dy od;
  a:=a+dx od;
end:

 
computelines:=proc(f,Ax,Ay,p,q,r,s,m,n)
local dx,dy,k,ldx,slp,ldy,xinc,yinc,a,b,i,j;
dx:=(q-p)/m; dy:=(s-r)/n;
ldx:=3*dx/8; ldy:=3*dy/8;
a:=p;
for i from 0 to m do
  b:=r;
  for j from 0 to n do
    k:=2*j;
    slp:=traperror(evalf(f(a,b)));
    if lasterror=slp then
      xinc:=0; yinc:=0;
    elif not type(slp,'numeric') then
      xinc:=0; yinc:=0;
    elif slp > dy/dx or -slp > dy/dx then
      xinc:=ldy/slp; yinc:=ldy;
    else xinc:=ldx; yinc:=ldx*slp fi;
    Ax[i,k]:=a-xinc;
    Ay[i,k]:=b-yinc;
    Ax[i,k+1]:=a+xinc;
    Ay[i,k+1]:=b+yinc;
    b:=b+dy od;
  a:=a+dx od;
end:

macro( INITCOND1 = {[constant,constant]} );
fieldplot:=proc(F,H,V)
local npts,init,flow,A,Ax,Ay,i,j,nsteps,f,numstyle,
 CC,h,m,n,p,q,r,s,cargs,flows,precision,step,
 fl, flines,carg,initlist,Aer:
options `Copyright 1991 by Daniel Schwalbe`;
precision := max(Digits,trunc(evalhf(Digits)));
npts:=1; m:=16; n:=12; nsteps:=20:
if not type(F,{procedure,numeric}) then
   ERROR(`1st argument (functions) should be a procedure`) fi;
if not type([H,V],[constant..constant,constant..constant]) then
   ERROR(`2nd and 3rd arguments (domains) must be real ranges`) fi;
step:= proc(x) if x<0 then 0 else 1 fi end;
if type(F,numeric) then f:=subs(_body=F,proc() _body end) else
  f:=subs(Heaviside=step,eval(F)) fi:
p:=evalf(lhs(H),precision); q:=evalf(rhs(H),precision);
r:=evalf(lhs(V),precision); s:=evalf(rhs(V),precision);
h:=evalf((q-p)/nsteps,precision):
initlist:={};
cargs:={args[4..nargs]};
for carg in cargs do
  if type(carg,'set') then
    if type(carg,'set(INITCOND1)') then initlist := carg
   else ERROR(`initial conditions must be a set of points of type`,INITCOND1)
    fi;
  elif type(carg,'list') then # A single initial condition
    if type(carg,INITCOND1) then initlist:={carg}
	else ERROR(`initial condition must be of type`,INITCOND1)
    fi;
  elif not type(carg,equation) then
	ERROR(`optional arguments must be equations`)
  elif op(1,carg)='stepsize' then 
     h:=evalf(op(2,carg),precision);
	 nsteps:=trunc((q-p)/h + 1);
     if not type(h,numeric) then
	ERROR(`stepsize must be a numerical value`)
     fi;
  elif op(1,carg)='numsteps' then
     nsteps:=op(2,carg);
     if not type(nsteps,posint) then
	ERROR(`numsteps must be a positive integer`)
     fi;
  elif op(1,carg)='iterations' then 
     npts:=op(2,carg);
     if not type(npts,posint) then
	ERROR(`number of iterations must be a positive integer`)
     fi;
  elif op(1,carg)='scheme' then numstyle := 'adaptive'
  elif op(1,carg)='grid' then
     if not type(rhs(carg),[integer,integer]) then
	ERROR(`grid must be a list of two integers`) fi;
     m:=op(2,carg)[1];
     n:=op(2,carg)[2]
  else ERROR(`unknown option`,lhs(carg))
  fi
od;
CC:=array(1..9);
CC[1]:=m: CC[2]:=n: CC[3]:=h:
CC[4]:=p: CC[5]:=q: CC[6]:=r: CC[7]:=s:
if m>0 and n>0 then
  Ax:=array(0..m,0..2*n+1);
  Ay:=array(0..m,0..2*n+1);
  fl:=traperror(evalhf(computelineshf(
            f,var(Ax),var(Ay),var(CC)))):
  if lasterror=fl then #print(`computelineshf failed`);
    computelines(f,Ax,Ay,p,q,r,s,m,n) fi;
  flines:=seq(seq(CURVE(BLACK,LINE,
      [Ax[i,2*j],Ay[i,2*j],Ax[i,2*j+1],
           Ay[i,2*j+1]]),i=0..m),j=0..n):
else flines:=[0,0] fi:
i:=1:
for init in initlist do
  CC[1]:=evalf(init[1],precision);
  CC[2]:=evalf(init[2],precision);
  if numstyle<>'adaptive' then
    A:=array(2*(-nsteps-3)..2*(nsteps)+3):
    Aer:=traperror(evalhf(fieldrk(var(A),f,var(CC),npts)));
    if Aer=lasterror then  #print(`evalhf(fieldrk) failed`);
      Aer:= traperror(fieldrk(A,f,CC,npts)) fi;
    if Aer<>lasterror then
      flow[i]:=seq(A[j],j=round(CC[8])..round(CC[9])):
      i:=i+1 fi;
   else
     A:=adaptiverk(f,CC):
     flow[i]:=seq(A[j],j=round(CC[8])..round(CC[9]));
     i:=i+1 fi od:

flows:=seq(CURVE(RED,LINE,[flow[j]]),j=1..i-1);
m:=max(m,16);n:=max(n,12);
PLOT(AXIS(HORIZONTAL,SCALE(4,LINEAR,NUMBER)),
 AXIS(VERTICAL,SCALE(4,LINEAR,NUMBER)),
 RANGE(HORIZONTAL,p-(q-p)/m..q+(q-p)/m),
 RANGE(VERTICAL,r-(s-r)/n..s+(s-r)/n),
 flines,flows):
end:

##################################################
rk2step:=proc(f1,f2,h,p,q,r,s,npts,pt) 
local t,x,y,k11,k12,k13,k14,k21,k22,k23,k24;
options `Copyright 1991 by Daniel Schwalbe`;
t:=pt[1]:x:=pt[2]:y:=pt[3]: 
to npts while x<=q and x>=p and y<=s and y>=r do
  k11:=f1(t,x,y):
  k21:=f2(t,x,y):
  k12:=f1(t+h/2,x+h*k11/2,y + h*k21/2):
  k22:=f2(t+h/2,x+h*k11/2,y + h*k21/2):
  k13:=f1(t+h/2,x+h*k12/2,y + h*k22/2):
  k23:=f2(t+h/2,x+h*k12/2,y + h*k22/2):
  k14:=f1(t+h,x+h*k13,y + h*k23):
  k24:=f2(t+h,x+h*k13,y + h*k23):
  x:=x+h*(k11+2*k12+2*k13+k14)/6;
  y:=y+h*(k21+2*k22+2*k23+k24)/6;
  t:=t+h; od;
pt[1]:=t: pt[2]:=x:pt[3]:=y:
end:

phaserk:=proc(A,f1,f2,CC,nsteps,npts)
local pta,ptb,pt,p,q,r,s,j,h;
options `Copyright 1991 by Daniel Schwalbe`;
pt:=array(1..3):
pta:=array(1..2):ptb:=array(1..2):
h:=CC[1]/npts:
pt[1]:=CC[2]: pt[2]:=CC[3]: pt[3]:=CC[4]:
A[0]:=pt[2] :A[1]:=pt[3]:
p:=CC[5]:q:=CC[6]:r:=CC[7]:s:=CC[8]:
j:=0:
while j<nsteps and pt[2]<=q and pt[2]>=p and 
  pt[3]<=s and pt[3]>=r do
    j:=j+1:
    rk2step(f1,f2,h,p,q,r,s,npts,pt):
    A[2*j]:=pt[2]: A[2*j+1]:=pt[3]: od:
if pt[2]>q or pt[2]<p or pt[3]>s or pt[3]<r then
  pta[1]:=A[2*j-2]:pta[2]:=A[2*j-1]:
  ptb[1]:=A[2*j]: ptb[2]:=A[2*j+1]:
  intbox(p,q,r,s,pta,ptb):
  A[2*j]:=ptb[1]:A[2*j+1]:=ptb[2] fi:
CC[10]:=2*j+1:
j:=0:
pt[1]:=CC[2]: pt[2]:=CC[3]: pt[3]:=CC[4]:
while -j<nsteps and pt[2]<=q and pt[2]>=p and
  pt[3]<=s and pt[3]>=r do
    j:=j-1:
    rk2step(f1,f2,-h,p,q,r,s,npts,pt):
    A[2*j]:=pt[2]: A[2*j+1]:=pt[3]: od:
if pt[2]>q or pt[2]<p or pt[3]>s or pt[3]<r then
  pta[1]:=A[2*j+2]:pta[2]:=A[2*j+3]:
  ptb[1]:=A[2*j]: ptb[2]:=A[2*j+1]:
  intbox(p,q,r,s,pta,ptb):
  A[2*j]:=ptb[1]:A[2*j+1]:=ptb[2] fi:
CC[9]:=2*j
end:


computevectshf:=proc(Ax,Ay,fa,fb,p,q,r,s,m,n)
local i,j,k,maxlen,len,a,b,x,y,ap,bp,
 scale,wx,wy,zx,zy,blocks,dx,dy;
options `Copyright 1991 by Daniel Schwalbe`;
  dx:=(q-p)/m: dy:=(s-r)/n:
if m > n then blocks:= m else blocks:= n fi;
maxlen:=0;
a:=p; 
for j from 0 to m do
  b:=r;
  for i from 0 to n do
    x:=fa(0,a,b);
    y:=fb(0,a,b);
    len:=((blocks/(q-p)*x)^2+(blocks/(s-r)*y)^2)^(1/2);
    if len > maxlen then maxlen:=len fi; 
      b:=b+dy od; a:=a+dx od;
a:=p;
if blocks < 20 then blocks:=20 fi;
for j from 0 to m do
  b:=r;
  for i from 0 to n do
    k:=6*i;
    x:=fa(0,a,b)/maxlen; y:=fb(0,a,b)/maxlen;
    if x=0 and y=0 then 
      scale:= 0 else 
      scale:=(((blocks*x/(q-p))^2+(blocks*y/(s-r))^2))^(-1/2) fi;
    wx:= scale*x;   wy:=scale*y; 
 #zx and zy are components of perpendicular vector 
    zx:= -scale*y/(s-r)*(q-p); zy:= scale*x/(q-p)*(s-r) ;
    Ax[j,k]:=a;
    Ay[j,k]:=b;
    ap:=a+2*x/3+wx/3; bp:=b+2*y/3+wy/3;
    Ax[j,k+1]:=ap; 
    Ay[j,k+1]:=bp;
    Ax[j,k+2]:=ap-wx/6+zx/6; 
    Ay[j,k+2]:=bp-wy/6+zy/6;
    Ax[j,k+3]:=ap;#+wx/6; 
    Ay[j,k+3]:=bp;#+wy/6;
    Ax[j,k+4]:=ap-wx/6-zx/6; 
    Ay[j,k+4]:=bp-wy/6-zy/6;
    Ax[j,k+5]:=ap; 
    Ay[j,k+5]:=bp;
    b:=b+dy od; a:=a+dx od
end:

#Only difference in computevects from computevectshf is trapping
#for errors when evaluating input procedures.

computevects:=proc(Ax,Ay,fa,fb,p,q,r,s,m,n)
local i,j,k,maxlen,len,a,b,x,y,ap,bp,
 scale,wx,wy,zx,zy,blocks,dx,dy;
options `Copyright 1991 by Daniel Schwalbe`;
  dx:=(q-p)/m: dy:=(s-r)/n:
if m > n then blocks:= m else blocks:= n fi;
maxlen:=0;
a:=p; 
for j from 0 to m do
  b:=r;
  for i from 0 to n do
    x:= traperror(fa(0,a,b));
	if x = lasterror then y:=0; x:=0 else
    y:=traperror(fb(0,a,b)) fi;
	if y = lasterror then y:=0; x:=0 fi;
	if not type([x,y],[numeric,numeric]) then
	  x:=0; y:=0 fi;
    len:=((blocks/(q-p)*x)^2+(blocks/(s-r)*y)^2)^(1/2);
    if len > maxlen then maxlen:=len fi; 
      b:=b+dy od; a:=a+dx od;
a:=p;
if blocks < 20 then blocks:=20 fi;
for j from 0 to m do
  b:=r;
  for i from 0 to n do
    k:=6*i;
    x:= traperror(fa(0,a,b)/maxlen);
	if x = lasterror then y:=0; x:=0 else
    y:=traperror(fb(0,a,b)/maxlen) fi;
	if y = lasterror then y:=0; x:=0 fi;
	if not type([x,y],[numeric,numeric]) then
	  x:=0; y:=0 fi;
    if x=0 and y=0 then 
      scale:= 0 else 
      scale:=(((blocks*x/(q-p))^2+(blocks*y/(s-r))^2))^(-1/2) fi;
    wx:= scale*x;   wy:=scale*y; 
 #zx and zy are components of perpendicular vector 
    zx:= -scale*y/(s-r)*(q-p); zy:= scale*x/(q-p)*(s-r) ;
    Ax[j,k]:=a;
    Ay[j,k]:=b;
    ap:=a+2*x/3+wx/3; bp:=b+2*y/3+wy/3;
    Ax[j,k+1]:=ap; 
    Ay[j,k+1]:=bp;
    Ax[j,k+2]:=ap-wx/6+zx/6; 
    Ay[j,k+2]:=bp-wy/6+zy/6;
    Ax[j,k+3]:=ap;#+wx/6; 
    Ay[j,k+3]:=bp;#+wy/6;
    Ax[j,k+4]:=ap-wx/6-zx/6; 
    Ay[j,k+4]:=bp-wy/6-zy/6;
    Ax[j,k+5]:=ap; 
    Ay[j,k+5]:=bp;
    b:=b+dy od; a:=a+dx od
end:



macro( INITCOND = {[constant,constant],[constant,constant,constant]} );
phaseplot:=proc(F,H,V)
local i,j,p,q,r,s,phasevects,initlist,precision,
  pv,init,A,flow,f1,f2,npts,nsteps,m,n,step,
  h,cargs,carg,Alist,Ax,Ay,CC,flows:
options `Copyright 1991 by Daniel Schwalbe`;
if not type(F,[{procedure,numeric},{procedure,numeric}]) then
   ERROR(`1st argument (functions) should be a list of two procedures`) fi;
if not type([H,V],[constant..constant,constant..constant]) then
   ERROR(`2nd and 3rd arguments (domains) must be real ranges`) fi;
step:=proc(x) if x<0 then 0 else 1 fi end:
nsteps:=50:npts:=1;h:=0.1:m:=16:n:=14:
precision := max(Digits,trunc(evalhf(Digits)));
if type(F[1],numeric) then f1:=subs(_body=F[1],proc() _body end) else
  f1:= subs(Heaviside=step,eval(F[1])) fi;
if type(F[2],numeric) then f2:=subs(_body=F[2],proc() _body end) else
  f2:= subs(Heaviside=step,eval(F[2])) fi;
p:=evalf(lhs(H),precision); q:=evalf(rhs(H),precision);
r:=evalf(lhs(V),precision); s:=evalf(rhs(V),precision);
cargs:={args[4..nargs]};
initlist:={};
for carg in cargs do
  if type(carg,'set') then
    if type(carg,'set(INITCOND)') then initlist := carg
    else ERROR(`initial conditions must be a set of points of type`,INITCOND)
    fi;
  elif type(carg,'list') then # A single initial condition
    if type(carg,INITCOND) then initlist:={carg}
	else ERROR(`initial condition must be of type`,INITCOND)
    fi;
  elif not type(carg,equation) then
	ERROR(`optional arguments must be equations`)
  elif op(1,carg)='stepsize' then 
     h:=evalf(op(2,carg),precision);
     if not type(h,numeric) then
	ERROR(`stepsize must be a numerical value`)
     fi;
  elif op(1,carg)='numsteps' then
     nsteps:=op(2,carg);
     if not type(nsteps,posint) then
	ERROR(`numsteps must be a positive integer`)
     fi;
  elif op(1,carg)='iterations' then 
     npts:=op(2,carg);
     if not type(npts,posint) then
	ERROR(`number of iterations must be a positive integer`)
     fi;
  elif op(1,carg)='grid' then
     if not type(rhs(carg),[integer,integer]) then
	ERROR(`grid must be a list of two positive integers`) fi;
     m:=op(2,carg)[1];
     n:=op(2,carg)[2]
  else ERROR(`unknown option`,lhs(carg))
  fi
od;

if m>0 and n>0 then
  Ax:=array(0..m,0..6*n+5);
  Ay:=array(0..m,0..6*n+5);
  pv:=traperror(evalhf(computevectshf(var(Ax),var(Ay),
     f1,f2,p,q,r,s,m,n)));
  if lasterror=pv then #print(`computevectshf failed`);
    computevects(Ax,Ay,f1,f2,p,q,r,s,m,n) fi;
  phasevects:=seq(seq(CURVE(BLACK,LINE,
      [Ax[j,6*i],Ay[j,6*i],Ax[j,6*i+1],Ay[j,6*i+1],
      Ax[j,6*i+2],Ay[j,6*i+2],Ax[j,6*i+3],Ay[j,6*i+3],
      Ax[j,6*i+4],Ay[j,6*i+4],Ax[j,6*i+5],Ay[j,6*i+5]]),
        j=0..m),i=0..n):
  else phasevects:=[0,0] fi;

i:=1;
for init in initlist do
  A:=array(-2*nsteps..2*nsteps+1);
  CC:=array(1..10):
  CC[1]:=h:
  if nops(init)=3 then
    CC[2]:=evalf(init[1]);
    CC[3]:=evalf(init[2]);
    CC[4]:=evalf(init[3]);
  else
    CC[2]:=0;
    CC[3]:=evalf(init[1]);
    CC[4]:=evalf(init[2]) fi;
  CC[5]:=p: CC[6]:=q: CC[7]:=r: CC[8]:=s:
  pv:=traperror(evalhf(phaserk(
  var(A),f1,f2,var(CC),nsteps,npts)));
  if pv = lasterror then
    phaserk(A,f1,f2,CC,nsteps,npts) fi;
  Alist:=seq(A[j],j=round(CC[9])..round(CC[10])):
  flow[i]:= CURVE(RED,LINE,[Alist]);
  i:=i+1; od; 

flows:=seq(flow[j],j=1..i-1):
m:=max(m,16);n:=max(n,12);
 
PLOT(AXIS(HORIZONTAL,SCALE(4,LINEAR,NUMBER)),
 AXIS(VERTICAL,SCALE(4,LINEAR,NUMBER)),
 RANGE(HORIZONTAL,p-(q-p)/m..q+(q-p)/m),
 RANGE(VERTICAL,r-(s-r)/n..s+(s-r)/n),
 phasevects,flows):
end:


##################################################

picard:=proc(f,init,n)
 local parray,ptemp,t0,y0:
options `Copyright 1991 by Daniel Schwalbe`;
 t0:=init[1]; y0:=init[2];
 ptemp:=y0+int(f('s',y0),'s'=t0..'t');
 parray:=ptemp;
 from 2 to n do
 ptemp:=y0+int(f('s',subs('t'='s',ptemp)),'s'=t0..'t');
 parray:=parray,sort(ptemp) od;
 parray;
end:

##################################################

#save `ODE.m`;
#quit
