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)