a = 0.5
b = 2
s = plot3d(x^2 + y^2, (x, -4, 4), (y, -4, 4), color='gray', opacity=0.5)
px = parametric_plot3d([x, b, z], (x, -4, 4), (z, 0, 9), color='red', opacity=0.5)
py = parametric_plot3d([a, y, z], (y, -4, 4), (z, 0, 9), color='blue', opacity=0.5)
cx = parametric_plot3d([x, b, x^2 + b^2], (x, -4, 4), color='black')
cy = parametric_plot3d([a, y, a^2 + y^2], (y, -4, 4), color='black')
dx = parametric_plot3d([x, b, 2*a*(x-a) + a^2 + b^2], (x, -4, 4), color='black')
dy = parametric_plot3d([a, y, 2*b*(y-b) + a^2 + b^2], (y, -4, 4), color='black')
show(s+px+py+cx+cy+dx+dy)
