I am bleary-eyed from trying to get Maple to follow simple syntactic transformations, so maybe I just don't see the obvious mistake in all this... I would expect I2 to equal I3 (below), since constants can be factored out of integrals: with(student); assume(m>0, mu>0, s>0, sp>0, y0 >0, y1>0, z0>0, z1>0, c0::real, c1::real); xi:=(y,z,e) -> exp(-z^2*(y-e)^2/(2*sp^2)); F:= (u,v) -> ((c0-c1*u)/s)*exp(-(u-mu)^2/(2*s^2) - (c0-c1*u)*v); Xi:= (u,v) -> xi(y0,z0,u)*xi(y1,z1,v); T:= (u,v) -> Xi(u,v)*F(u,v); alpha:= z0^2/(2*sp^2); beta:= z1^2/(2*sp^2); gama:= 1/(2*s^2); epsilon:= u*c1/(2*beta) + (2*beta*y1-c0)/(2*beta); delta:= z1/(2^(1/2)*sp); Fu:= exp(-alpha*(y0-u)^2 - gama*(u-mu)^2 + beta*epsilon^2 - beta*y1^2)*((c0-c1*u)/s); Fuv:= exp(-beta*(v-epsilon)^2); simplify(T(u,v) - Fu*Fuv); I0:= int(Fu*(1/delta)*int(exp(-x^2),x=-delta*epsilon..infinity),u=0..m); simplify(value(I0 - Doubleint(T(u,v),v=0..infinity,u=0..m))); phi:= 1/2*z0^2/sp^2+1/2/s^2-1/2/z1^2*sp^2*c1^2; theta:= (c1*y1+1/s^2*mu-sp^2/z1^2*c1*c0+z0^2/sp^2*y0)/(2*phi); psi:= -1/2*z0^2/sp^2*y0^2-1/2/s^2*mu^2-y1*c0+1/2*sp^2/z1^2*c0^2; I1:= int(exp(-phi*(u-theta)^2+phi*theta^2+psi)*((c0-c1*u)/s)*int(exp(-x^2),x=-delta*epsilon..infinity),u=0..m)/delta; simplify(value(I1-I0)); I2:= int(exp(phi*theta^2+psi)*exp(-phi*(u-theta)^2)*((c0-c1*u)/s)*int(exp(-x^2),x=-delta*epsilon..infinity),u=0..m)/delta; simplify(value(I2-I1)); diff(exp(phi*theta^2+psi),u); I3:= exp(phi*theta^2+psi)*int(exp(-phi*(u-theta)^2)*((c0-c1*u)/s)*int(exp(-x^2),x=-delta*epsilon..infinity),u=0..m)/delta; simplify(value(I3-I2));

Appreciate the assume magic, but I don't have your mojo... with(student); assume(m>0, mu>0, s>0, sp>0, y0::real, y1::real, z0::real, z1::real, c0::real, c1::real); xi:=(y,z,e) -> exp(-z^2*(y-e)^2/(2*sp^2)); F:= (u,v) -> ((c0-c1*u)/s)*exp(-(u-mu)^2/(2*s^2) - (c0-c1*u)*v); Xi:= (u,v) -> xi(y0,z0,u)*xi(y1,z1,v); T:= (u,v) -> Xi(u,v)*F(u,v); simplify(int(int(T(u,v),v=0..infinity),u=0..m) - Doubleint(T(u,v),v=0..infinity,u=0..m)); simplify(int(int(q(u,v),v=0..infinity),u=0..m) - Doubleint(q(u,v),v=0..infinity,u=0..m));

alpha:= z0^2/(2*sp^2); beta:= z1^2/(2*sp^2); gama:= 1/(2*s^2); epsilon:= u*c1/(2*beta) + (2*beta*y1-c0)/(2*beta); delta:= z1/(2^(1/2)*sp); Fu:= exp(-alpha*(y0-u)^2 - gama*(u-mu)^2 + beta*epsilon^2 - beta*y1^2)*((c0-c1*u)/s); Fuv:= exp(-beta*(v-epsilon)^2); I1:= int(int(Fu*Fuv,v=0..infinity),u=0..max); I2:= int(Fu*int(Fuv,v=0..infinity),u=0..max); simplify(I2-I1);