The problem lies in the Maxima to Sage backtranslation :
sage: var('x','y','x_o','y_o','a','b', domain="real")
(x, y, x_o, y_o, a, b)
sage: eq02=y - y_o == -b^2*(x - x_o)*x/(a^2*y)
sage: eq2=x^2/a^2 + y^2/b^2 == 1
sage: solve([eq2,eq02],[x,y])
[ Snip.. ]
TypeError: unable to make sense of Maxima expression '[if((_SAGE_VAR_a^2*_SAGE_VAR_b^2!=0)and(-(_SAGE_VAR_a^2*(_SAGE_VAR_b^2*_SAGE_VAR_x_o*sqrt(_SAGE_VAR_a^2*_SAGE_VAR_y_o^2+_SAGE_VAR_b^2*_SAGE_VAR_x_o^2-_SAGE_VAR_a^2*_SAGE_VAR_b^2)-_SAGE_VAR_a^2*_SAGE_VAR_b^2*_SAGE_VAR_y_o))/(_SAGE_VAR_a^2*_SAGE_VAR_y_o^2+_SAGE_VAR_b^2*_SAGE_VAR_x_o^2)!=0),[_SAGE_VAR_x==(_SAGE_VAR_a^2*_SAGE_VAR_y_o*sqrt(_SAGE_VAR_a^2*_SAGE_VAR_y_o^2+_SAGE_VAR_b^2*_SAGE_VAR_x_o^2-_SAGE_VAR_a^2*_SAGE_VAR_b^2)+_SAGE_VAR_a^2*_SAGE_VAR_b^2*_SAGE_VAR_x_o)/(_SAGE_VAR_a^2*_SAGE_VAR_y_o^2+_SAGE_VAR_b^2*_SAGE_VAR_x_o^2),_SAGE_VAR_y==(_SAGE_VAR_a^2*_SAGE_VAR_b^2*_SAGE_VAR_y_o-_SAGE_VAR_b^2*_SAGE_VAR_x_o*sqrt(_SAGE_VAR_a^2*_SAGE_VAR_y_o^2+_SAGE_VAR_b^2*_SAGE_VAR_x_o^2-_SAGE_VAR_a^2*_SAGE_VAR_b^2))/(_SAGE_VAR_a^2*_SAGE_VAR_y_o^2+_SAGE_VAR_b^2*_SAGE_VAR_x_o^2)],union()),if((_SAGE_VAR_a^2*_SAGE_VAR_b^2!=0)and((_SAGE_VAR_a^2*(_SAGE_VAR_b^2*_SAGE_VAR_x_o*sqrt(_SAGE_VAR_a^2*_SAGE_VAR_y_o^2+_SAGE_VAR_b^2*_SAGE_VAR_x_o^2-_SAGE_VAR_a^2*_SAGE_VAR_b^2)+_SAGE_VAR_a^2*_SAGE_VAR_b^2*_SAGE_VAR_y_o))/(_SAGE_VAR_a^2*_SAGE_VAR_y_o^2+_SAGE_VAR_b^2*_SAGE_VAR_x_o^2)!=0),[_SAGE_VAR_x==(_SAGE_VAR_a^2*_SAGE_VAR_b^2*_SAGE_VAR_x_o-_SAGE_VAR_a^2*_SAGE_VAR_y_o*sqrt(_SAGE_VAR_a^2*_SAGE_VAR_y_o^2+_SAGE_VAR_b^2*_SAGE_VAR_x_o^2-_SAGE_VAR_a^2*_SAGE_VAR_b^2))/(_SAGE_VAR_a^2*_SAGE_VAR_y_o^2+_SAGE_VAR_b^2*_SAGE_VAR_x_o^2),_SAGE_VAR_y==(_SAGE_VAR_b^2*_SAGE_VAR_x_o*sqrt(_SAGE_VAR_a^2*_SAGE_VAR_y_o^2+_SAGE_VAR_b^2*_SAGE_VAR_x_o^2-_SAGE_VAR_a^2*_SAGE_VAR_b^2)+_SAGE_VAR_a^2*_SAGE_VAR_b^2*_SAGE_VAR_y_o)/(_SAGE_VAR_a^2*_SAGE_VAR_y_o^2+_SAGE_VAR_b^2*_SAGE_VAR_x_o^2)],union())]' in Sage
and indeed :
sage: SolM=solve([eq2,(eq02*y/b^2).expand()],[x,y],solution_dict=True) ; SolM
[{x: (a^2*b^2*x_o - sqrt(-a^2*b^2 + b^2*x_o^2 + a^2*y_o^2)*a^2*y_o)/(b^2*x_o^2 + a^2*y_o^2),
y: (a^2*b^2*y_o + sqrt(-a^2*b^2 + b^2*x_o^2 + a^2*y_o^2)*b^2*x_o)/(b^2*x_o^2 + a^2*y_o^2)},
{x: (a^2*b^2*x_o + sqrt(-a^2*b^2 + b^2*x_o^2 + a^2*y_o^2)*a^2*y_o)/(b^2*x_o^2 + a^2*y_o^2),
y: (a^2*b^2*y_o - sqrt(-a^2*b^2 + b^2*x_o^2 + a^2*y_o^2)*b^2*x_o)/(b^2*x_o^2 + a^2*y_o^2)}]
sage: SolP=solve([eq2,(eq02*y/b^2).expand()],[x,y],algorithme="sympy",solution_dict=True) ; SolP
[{x: (a^2*b^2*x_o - sqrt(-a^2*b^2 + b^2*x_o^2 + a^2*y_o^2)*a^2*y_o)/(b^2*x_o^2 + a^2*y_o^2),
y: (a^2*b^2*y_o + sqrt(-a^2*b^2 + b^2*x_o^2 + a^2*y_o^2)*b^2*x_o)/(b^2*x_o^2 + a^2*y_o^2)},
{x: (a^2*b^2*x_o + sqrt(-a^2*b^2 + b^2*x_o^2 + a^2*y_o^2)*a^2*y_o)/(b^2*x_o^2 + a^2*y_o^2),
y: (a^2*b^2*y_o - sqrt(-a^2*b^2 + b^2*x_o^2 + a^2*y_o^2)*b^2*x_o)/(b^2*x_o^2 + a^2*y_o^2)}]
sage: bool(Set(SolM)==Set(SolP))
True
That's probably Trac#13773.
Shit happens. Complain fiercely...
HTH,