def f(x): return x/(x**2+1) ymin=f(x).subs(x,const[0]) ymax=f(x).subs(x,const[1]) from scipy.optimize import fmin xmin=fmin(f,0) xmax=fmin(lambda x: -f(x),0) integ=integrate(f(x),(x,xmin,xmax)) u,v=symbols("u,v") x=cos(u)*(cos(v)+3) y=sin(u)*(cos(v)+3) z=sin(v)+u plot3d_parametric_surface(x,y,z,(u,umin,umax),(v,vmin,vmax),...) solut=dsolve(difeq,T(x)) const=solve([solut.rhs.subs(x,0)-T1,solut.rhs.subs(x,s)-T2]) funct=solut.subs(const)