var('x', 'y', 'z')
f(x, y) = x^2 + y^2
a = 2
b = 1
fx(x, y) = f(x, y).derivative(x)
fy(x, y) = f(x, y).derivative(y)
S = plot3d(f(x,y), (x, -3, 3), (y, -3, 3), color = 'blue', opacity = 0.5)
Px = parametric_plot3d([x, b, z], (x, -3, 3), (z, 0, 10), color = 'red', opacity = 0.5)
Cx = parametric_plot3d([x, b, f(x, b)], (x, -3, 3), color = 'black')
Tx = parametric_plot3d([a + x, b, f(a, b) + fx(a, b)*x], (x, -3 - a, 3 - a), color = 'black')
Py = parametric_plot3d([a, y, z], (y, -3, 3), (z, 0, 10), color = 'magenta', opacity = 0.5)
Cy = parametric_plot3d([a, y, f(a, y)], (y, -3, 3), color = 'black')
Ty = parametric_plot3d([a, b + y, f(a, b) + fy(a, b)*y], (y, -3 - b, 3 - b), color = 'black')
P = plot3d(f(a,b) + fx(a, b)*(x-a) + fy(a, b)*(y-b), (x, -3, 3), (y, -3, 3), color = 'green')
show(S + Ty + Cy + Py, aspect_ratio = [1, 1, 0.5])
