Estoy intentando hacer una integral doble en fortran por el método de simpson pero no se como declarar los limites de integracion.El problema dice asi:
Calcular mediante la formula del Simpson la
integral de x*exp(-x*y) siendo la
region (0 ≤ x ≤ 2, u(x) ≤ y ≤ v(x)} y siendo
u(x) = x**2-2*x+2 y
v(x) =3*exp(x) (para cada valor de x, c = u(x) y d = v(x)). Tomar n = 200 y utilizar dos subprogramas
funcion para calcular u(x) y v(x).
yo he hecho esto y me da error pero no se donde esta
program prosim
implicit none
integer::n,i
real*8::x,a,b,c,d,s,s1,s2,h
real*8,external::f,g_simp,u,v
n=256
a=0d0
b=2d0
c=u(x)
d=v(x)
h=(b-a)/dble(n)
s1=0d0
do i=1,n-1,2
x=a+dble(i)*h
s1=s1+g_simp(x,c,d,h)
end do
s2=0d0
do i=2,n-2,2
x=a+dble(i)*h
s2=s2+g_simp(x,c,d,h)
end do
s=(h/3d0)*(g_simp(a,c,d,h)+g_simp(b,c,d,h)+4*s1+2*s2)
print*, s
end program prosim
function g_simp(x,c,d,h)
implicit none
integer::j,m
real*8::h,s1,s2,c,d,k,g_simp,y,x
real*8,external::f,
m=2*int((d-c+2*h)/(2*h))
k=(d-c)/dble(m)
s1=0d0
do j=1,m-1,2
y=c+dble(j)*k
s1=s1+f(x,y)
end do
s2=0d0
do j=2,m-2,2
y=c+dble(j)*k
s2=s2+f(x,y)
end do
g_simp=(k/3d0)*(f(x,c)+f(x,d)+4*s1+2*s2)
end function g_simp
function f(x,y)
real*8::f,x,y
f=x*exp((-x)*y)
end function
function u(x)
real*8::u,x
u=x**2-2*x+2
end function
function v(x)
real*8::v,x
v=3*exp(x)
end function
