eqn =

NO. You think there is a unique solution. But that is wrong.

syms a1 a2 a3 a4 x; % Define symbolic variables

% Define the equation with all terms on one side (zero on the other side)

eqn = a1 + a4*(x^3 + x^2 + x + 1) + a3*(x^2 + x + 1) + a2*(x + 1) - x^2 == 0;

eqn = expand(eqn)

For that equation to hold true, for ALL values of x, this is what solve returns. What makes that equation valid for EVERY POSSIBLE value of x? That is, identically EQUAL TO ZERO.

[A, B] = equationsToMatrix(eqn, [a1, a2, a3, a4])

% Solve the linear system

sol = linsolve(A, B);

% Display the solution

disp('The solutions for a1, a2, a3, a4 are:');

disp(sol);

It holds true for ALL x, only when a1 = x^2, and all of the other coefficients are zero.

In fact, linsolve tells us there are other solutions that are possible. That is, the solution was not a unique choice. There is in fact a 3 dimensional solution space of infinitely many solutions.

solgen = solve(eqn,{a1,a2,a3,a4},returnconditions = true)

Solve tells us that, although it is not very clear what that space looks like. Effectively, it is controlled by the equation itself. So you could choose ANY value for z there, and z1, and z2.

syms z z1 z2

sol124 = subs(solgen,[z,z1,z2],[1 2 4])

Now we can plug that back into eqn...

subs(eqn,sol124)

and we see it is happy.

Now, COULD the solution you think is the correct one also work? Yes. But why do you think it is the ONLY unique solution? That is flat out wrong.

subs(eqn,[a1,a2,a3,a4],[0 -1 1 0])

Again, it is happy. BUT the solution linsolve returns is as equally valid.

Hmm, I think maybe what you want, is to see if a solution exists where all of the parameters a1,a2,a3,a4 are constant, and not a function of x. linsolve does not understand that as a requirement. Try this instead:

expr = a1 + a4*(x^3 + x^2 + x + 1) + a3*(x^2 + x + 1) + a2*(x + 1) - x^2

So I wrote it as an expression, instead of an equation. That allows me to use coeffs.

coeffs(expr,x) == 0

do you see that generates a set of equations for the coefficients of powers of x. Now we can use solve.

solve(coeffs(expr,x) == 0)