I get the correct result (at least for coe = [1;1;1;1;1;1]) :
aex = expandA(ones(6,1)).'
function aex = expandA(coe)
A = coe(1);
B = coe(2);
C = coe(3);
D = coe(4);
E = coe(5);
F = coe(6);
a_scale = 1.3;
b_scale = 1.3;
x0 = (2*C*D - B*E) / (B^2 - 4*A*C);
y0 = (2*A*E - B*D) / (B^2 - 4*A*C);
theta = 0.5 * atan2(-B, C - A);
numerator1 = 2*(A*E^2+C*D^2-B*D*E+(B^2-4*A*C)*F);
numerator2_1 = A+C+((A-C)^2+B^2)^0.5;
numerator2_2 = A+C-((A-C)^2+B^2)^0.5;
denominator = B^2-4*A*C;
a = -sqrt(numerator1*numerator2_1)/ denominator;
b = -sqrt(numerator1*numerator2_2)/ denominator;
ae = a_scale*a;
be = b_scale*b;
syms x y
x_prime = (x - x0) * cos(theta) + (y - y0) * sin(theta);
y_prime = -(x - x0) * sin(theta) + (y - y0) * cos(theta);
new_ellipse_eq = (be^2*x_prime^2) + (ae^2*y_prime^2)-ae^2*be^2;
new_ellipse_general = expand(new_ellipse_eq);
[coeffsd,t] = coeffs(new_ellipse_general, [x y]);
double(coeffsd)
% aex=zeros(6,1);
% aex(1)=double(coeffsd(1));
% aex(2)=double(coeffsd(2));
% aex(3)=double(coeffsd(4));
% aex(4)=double(coeffsd(3));
% aex(5)=double(coeffsd(5));
% aex(6)=double(coeffsd(6));
A2 = ae^2*sin(theta)^2+be^2*cos(theta)^2;
B2 = 2*(be^2-ae^2)*sin(theta)*cos(theta);
C2 = ae^2*cos(theta)^2+be^2*sin(theta)^2;
D2 = -2*A2*x0-B2*y0;
E2 = -B2*x0-2*C2*y0;
F2 = A2*x0^2+B2*x0*y0+C2*y0^2-ae^2*be^2;
aex = [A2;B2;C2;D2;E2;F2];
end
Which MATLAB version do you use ?
