2021-02-15 15:56:03 +0100 commented answer Multivariable numerical-symbolic integration Nope, still not working. 2021-02-14 18:29:53 +0100 commented question Multivariable numerical-symbolic integration From what I've understood on the docs numerical_integral forces you to set the value of every other non integrating variable in your integrand. I don't need a number, I possibly need an expression. 2021-02-14 18:04:54 +0100 commented answer Multivariable numerical-symbolic integration I didn't know the integration variables order mattered in sage, thanks! My integrand complies with Fubini's theorem. I'm giving your suggestion a try now but it's still really slow and it could solve my use case but not similar ones in the future. If it won't work I'll try to simplify my integral by substitution to avoid composite functions in the denominator which is the reason why the number of terms got multiplied by the number of points after every integration with earlier method. 2021-02-13 16:22:32 +0100 received badge ● Editor (source) 2021-02-13 16:22:04 +0100 commented question Multivariable numerical-symbolic integration No, I didn't mean that. My bad though, I've not explained clearly enough, I just edited the example code to make it clearer and added a (not working in my case) not built in solution to better explain what I mean. 2021-02-13 11:43:08 +0100 asked a question Multivariable numerical-symbolic integration I need to find the symbolic expression of a multivariable integral. The real integrand contains 11 variables (I have to integrate "only" four of them) so example code is simplified here: x, y, z = var('x','y','z') integrand = sin(x)*cos(y)*z result = integrand.integral(x, 0, pi).integral(z, 0, 1)  It works just fine generally but given the complexity of the integral no algorithm gives me a result and it returns just: integrate(integrate(sin(x)*cos(y), x, 0, pi), z, 0, 1)  Is there a built in way to numerically evaluate an integral on some variables that returns symbolic result with the numerical integration found coefficients (or constants)? Namely, in this case: cos(y)  EDIT1: I managed to get a very very simple example of the function I need: def numsym_integral( f, variable, a, b, points = 1000): dx = (b-a)/points result = 0 for i in range(0,points): result = result + f.substitute({variable:a + i*dx}) * dx return result  It still doesn't work though in my specific case because of the complexity of the expression. I get the number of terms gets multiplied by "points" per integration, leading to something too big. Could something built in be more efficient?