var('z,theta,r')
f(x,y,z)=(x,y,z*(9-x^2-y^2))
solid_opts=dict(color='yellow',opacity=0.3)
solid_graph=parametric_plot3d((r*cos(theta),r*sin(theta),9-r^2), (r,0,3),(theta,-pi/2,pi/2),**solid_opts)
#solid_graph+=parametric_plot((0,y,z*(9-y^2)), (y,-3,3),(z,0,1),**solid_opts)
solid_graph+=parametric_plot((x,y*sqrt(9-x^2),0), (x,0,3),(y,-1,1),**solid_opts)
cross_section_graph=sum(parametric_plot(f(x,i,z), (x,0,sqrt(9-i^2)),(z,0,1)) for i in [-2..2])
show(solid_graph+cross_section_graph)
