# computes error function (see Abramowitz/Stegun 7.1.26)
function erf(x)
dimension a(5)
data p/.3275911/
data a/.25482 9592,
      -.28449 6736,
      1.42141 3741,
     -1.45315 2027,
      1.06140 5429/

t = 1./(1+p*x)
erf = 0.
z = t
do i=1,5
{
  erf = erf + a(i)*z
  z = z*t
}
return(1 - exp(-x*x)*erf)
end
# main program
# computes Y' OF Y=erf*erc(xi*x) for x=[xl,xup]
write (1,104)
104 format('  eps')
read (1,100) eps  
100 format(f10.6)
write (1,101) 
101 format( ' xlow ')
read (1,100) xl
write (1,102)
102 format(' xup   ')
read (1,100) xu
repeat [
write (1,105)
105 format( ' xilow ')
read (1,100) xil
write (1,106)
106 format(' xiup   ')
read (1,100) xiu
write (1,107)
107 format(' xistep ')
read (1,100) xistep
xpage = 480/(xiu-xil)
ypage = 220/2.
call clgr
call initgr
ipen = 1
for (xi=xil;xi<=xiu;xi=xi+xistep)
       {
	x = zero(xl,xu,eps,xi)
#	write (1,109) xi,x
#	109 format('  xi ',f10.4,' Z ',f10.4)
	ixpage = ifix((xi-xil)*xpage + 30)
	iypage = ifix(x*ypage + 30)
	call plot (ixpage,iypage,ipen)
	ipen = 0
	xu = x
	}
]
end


function zero(xlow,xupp,eps,xi)
# computes the zero of y' using bisection
# needs the interval [xl,xu] and xi, returns
# z, where abs(y'(z)) <= eps
# returns -1 if xl and xu have y' of the same sign
#
	xl = xlow; xu = xupp
	if (der(xl,xi)*der(xu,xi) >> 0) return (-1)
	while (abs(xl-xu)>=eps)
	  {
	  x = (xl+xu)/2
	  z = der(x,xi)
	  if (z<<0) 
		xu = x 
	  else
		xl = x	
	  }
	return (x)
end
# der returns derivative of the dfunctional
# given x,xi
	function der(x,xi)
	  x1 = x	 
	  x2 = xi*x
	  z1 = erf(x1)
	  z2 = erf(x2)
	  z = (1-erf(x2))*exp(-x1**2)
	  z = z - xi*erf(x1)*exp(-x2**2)
	  return (z)
	end
