n1 = vector([1,-1,2]) #A normal vector for the first plane

p1 = vector([1,1,1]) #A point on the first plane

n2 = vector([2,1,1]) #A normal vector for the second plane

p2 = vector([-2,1,0]) #A point on the second plane

var('x y z')

r = vector([x,y,z])

eq1= n1.dot_product(r-p1)==0 #Equation of the first plane

eq2= n2.dot_product(r-p2)==0 #Equation of the second plane

pl1 = implicit_plot3d(eq1, (x,-5,5), (y,-5,5), (z,-5,5), color='red')

pl2 = implicit_plot3d(eq2, (x,-5,5), (y,-5,5), (z,-5,5), color='blue')

sol = solve([eq1,eq2], x, y)[0]

vecsol = vector([sol[0].rhs(), sol[1].rhs(), z]) #The line of intersection, using z as parameter

inter = parametric_plot3d(vecsol, (z,-5,5), color='magenta', thickness=8)

show(pl1+pl2+inter, axes=True, aspect_ratio=1)