c Homemade routine to integrate the integral equation:
c phi(x) = (sin(x)/2) + (sin(x)/pi)*Integral[dt sin(t) phi(t)] 
	parameter(N=1000)
	real phi(N),phi1(N),x(N)

c define pi=3.14..., initialize x(i):

	pi=acos(-1.0)
	do i=1,N
	x(i)=pi*(float(i)/float(N))
	phi(i)=1.0 
	end do

c Neumann iteration loop:
	do ii=1,3

c loop over x(i):
	do i=1,N

c loop over x'(i)"
	sum=0.0
	do j=1,N
	sum=sum+pi*sin(x(j))*phi(j)/float(N)
	end do
	phi1(i)=0.5*sin(x(i)) + sin(x(i))*sum/pi
	end do
	
c update phi(x):
	do i=1,N
	phi(i)=phi1(i)
	end do

	do i=1,N
	print*,x(i),phi1(i)
	end do
	print*,' '

	end do


	stop
	end