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/2,z), (x,0,sqrt(9-(i/2)^2)),(z,0,1)) for i in [-5..5])
cross_section_graph+=sum(parametric_plot(vector(f(x,i/2,z))+vector((0,1/2,0)), (x,0,sqrt(9-(i/2)^2)),(z,0,1)) for i in [-5..5])
cross_section_graph+=sum(parametric_plot(vector(f(x,i/2,1))+vector((0,y,0)), (x,0,sqrt(9-(i/2)^2)),(y,0,1/2)) for i in [-5..5])
show(solid_graph+cross_section_graph)
