@_Maxim_ 

I think you forgot index=1 in RootOf. Then it works as expected.
(It's only evalf which takes index=1 by default.)

@Mariusz Iwaniuk 

You cannot obtain this conclusion! The series could be divergent. for some c's.
As it is easy to check, the equation has at most 3 real roots,

@asa12 

Consider the row operation Eij of adding the row i to row j.
If  Eij . X = X  for each i,j (i<>j) then obviously X = O.

@asa12 

Your definition of "invariant matrix" is missing. For the standard definition, the only matrix is zero.

@_Maxim_ 

But if f is just absolutely integrable (in the Lebesgue sense) then the series could be divergent at any point! And anyway the side limits f(x0+0), f(x0-0) may not exist.

@sand15 

It works for me.

Have you loaded plots?
Or, use plots:-display

@tomleslie 

For Digits:=21 the imaginary part is still present.

@asa12 

You have successfully found what? The zero matrix or something else?

@asa12 

If you are considering row additions E, then a matrix X is invariant under them (i.e. E.X = X, for each E)
iff X = 0. Probably this is not what you need.

@asa12 

You should explain what you are trying to do. Code only (which does not work correctly) does not say much.

Are you looking for these?

https://en.wikipedia.org/wiki/Elementary_matrix

@_Maxim_ 

Actually, if f and f' are continuous except for a finite number of finite jumps then f has bounded variation and so the Fourier series is always convergent at x to the midpoint (f(x+)+f(x-))/2.

tsunamiBTP
for the rectangle in the open interval (-1,1)
f:=piecewise(x < -1/2, 0, x < 1/2, 1, 1/2 < x , 0);
the sum of the series (n=infinity) is
F:=piecewise(x < -1/2, 0, x=-1/2, 1/2, x < 1/2, 1, x=1/2, 1/2, 1/2 < x , 0):

I have just explained that.

@Markiyan Hirnyk 

Actually the infinity of solutions refers to ode=0, where:

ode:=(1/2)*(diff(y(x), x))^2 - (1-ln(y(x)^2))*y(x)^2;

This should be normal because the  solve(...)  was used and it gives 2 solutions.

Ya verifies this ode:

simplify( eval(ode, y(x)=Ya) );
    0

If we restrict to solve()[1]  then we have only two solutions.