Ok no worries thanks anyway.

I'll give it one more explanation and if still no replies I guess I'm doomed to solve this alone...

Ignore everything I've said and just concentrate on this.

How do you solve this equation?

(a + b)x + (a - b) = 3x + 4?

You compare coefficients!

We have,

(a + b) = 3 and

(a - b) = 4

This should give a = 3.5 and b = 0.5.

this is basically what I'm trying to do. Compare coefficients from two functions to find the value of the constants.

In my case, the functions are,

f(x) = 1 - lambda*x^2

and

alpha*f(f(x/alpha))

where alpha = 1/f(1) = 1/(1 - lambda)

Now the thing is, when you substitute f(x) = 1 - lambda*x^2 into alpha*f(f(x/alpha)) you get a function which also has x^4 terms which we do not have in f(x). So we just discard and ignore these higher order terms.

So working with the functions as above, we get,

alpha*f(f(x/alpha)) = 1 + (2*lambda^2 - 2*lambda^3)*x^2 + terms of order x^4 which we ignore.

So now we compare this to f(x) = 1 - lambda*x^2.

We then see that in order for the two functions to be equal (again, ignoring the x^4 terms), we must solve,

(2*lambda^2 - 2*lambda^3) = - lambda

This gives us lambda = (1 + sqrt(3))/2

We then substiture this back into out definition of alpha to get alpha being roughly -2.7.

This is actually a known constant (one of Feigenbaums constants) for which the true value is -2.5... so we are close but can improve.

Now! We add on extra terms to f(x) and repeat the above proceedure except now we will have more equations to solve with more unknowns.

So if we let f(x) = 1 - lambda*x^2 + a*x^4.

We now have alpha = 1/f(1) = 1/(1 - lambda + a)

Now if we were to write out alpha*f(f(x/alpha) in full, we'd see we have terms up x^16. We then discard anything above our highest term in f(x) which is x^4.

Now we will have coefficients in alpha*f(f(x/alpha)) we can compare to coefficients in f(x) and solve for a and lambda and thus find a value for alpha.

So! That's what I'm doing.

The more terms you add onto f(x), the more accurate alpha will get.

This is the code I'm currently using

with(RealDomain):

Order := 12:

f := x -> 1 - lambda*x^2 + a*x^4 + b*x^6 + c*x^8 + e*x^10:

alpha := 1/(1 - lambda + a + b + c + e):

Lambda := simplify(coeff(series(alpha*f(f(x/alpha)),x),x^2)):

A := simplify(coeff(series(alpha*f(f(x/alpha)),x), x^4)):

B := simplify(coeff(series(alpha*f(f(x/alpha)),x), x^6)):

C := simplify(coeff(series(alpha*f(f(x/alpha)),x), x^8)):

E := simplify(coeff(series(alpha*f(f(x/alpha)),x), x^10)):

Eq := {Lambda = -lambda,A = a,B = b,C = c,E = e}:

fsolve(Eq,{lambda,a,b,c,e},avoid={lambda=0,a=0,b=0,c=0,e=0});

But it will only go as high as x*10. Have left on a version where I had x^12 as the highest term but never computed even after a day of running.

So, if you understood this... Any idea how I can make my code more efficient.