# 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 erf*erc(xi*x) for x=[xl,xup]
repeat {
	write (1,101) 
	101 format( ' xlow ')
	read (1,100) xl
	write (1,102)
	102 format(' xup   ')
	read (1,100) xup
	write (1,103)
	103 format(' ksi ')
	read (1,100) xi
	for (x=xl; x<=xup; x=x+.01)
	  {
	  Z = ERF(X)
	  z = z*(1-erf(xi*x))
	  write (1,100) x,z
	  100 format(2f10.6)
	  }
}
end
