with(geometry):

nestedPolygons := (n)->draw([seq(rotation(poly||i, RegularPolygon(``, 2^(i+1), point(``, [0,0]), 1/cos(Pi/2^(i+1))), Pi/2^(i+1), `counterclockwise`), i=1..n)], axes=`NONE`):

nestedPolygons(5);

Not sure whether I've understood your problem correctly, but if xi represents a real (or complex) number then try:

w := (x)->a[0]+sum(a[i]*x^i, i=1..4);     # define the polynom

solve([w(0), w(1), D(D(w))(0), D(D(w))(1)]);   # solve system of equations, D() is the differential operator