program ily
real::a,b,n,h,p,i

open(7,file='int.dat',status='new')
a=-0.4
b=0.4
n=1000
p=a
h=(b-a)/n
i=(h/2)*(f(a)+f(b))

1 p=p+h
i=i+(h/2)*2*f(p+h)
write(7,*) p,i
if (p<=b) goto 1

write(*,*) i

end program ily

function f(x)
real::r,pi
r=0.1 ; pi=acos(-1.)
f=(6-5*x**5)*(exp(-(x**2/r**2))/(r*sqrt(pi)))
end function f