**Partial (non-)answer,** to be completed if divine inspiration strikes ;-)...

Solving this system in `SR`

doesn't seem to be possible.

```
sage: var("a, b, c, d")
(a, b, c, d)
sage: eq1 = c==a*(a/(a+2*b))+2*b*(a/(a+b)) ; eq2 = d==a*(b/(a+2*b))+b*(b/(a+b))
```

Solving the first equation for `b`

is trivial, leading to two solutions

```
sage: Solb1 = solve(eq1, b, solution_dict=True)
sage: Solb1
[{b: -1/4*(3*a^2 - 3*a*c + sqrt(-7*a^2 + 6*a*c + c^2)*a)/(2*a - c)},
{b: -1/4*(3*a^2 - 3*a*c - sqrt(-7*a^2 + 6*a*c + c^2)*a)/(2*a - c)}]
```

Substituting this solution in the second equation leads to non-linear equations, e. g. :

```
sage: eq3 = eq2.subs(Solb1[0])
sage: eq3
d == -1/2*(3*a^2 - 3*a*c + sqrt(-7*a^2 + 6*a*c + c^2)*a)*a/((2*a - c)*(2*a - (3*a^2 - 3*a*c + sqrt(-7*a^2 + 6*a*c + c^2)*a)/(2*a - c))) + 1/4*(3*a^2 - 3*a*c + sqrt(-7*a^2 + 6*a*c + c^2)*a)^2/((4*a - (3*a^2 - 3*a*c + sqrt(-7*a^2 + 6*a*c + c^2)*a)/(2*a - c))*(2*a - c)^2)
```

This equation cannot be solved by Maxima's solver : calling it without extra arguments fails :

```
sage: Solam = solve(eq3, a, algorithm="maxima", solution_dict=True) # Apparent success...
sage: len(Solam)
4
```

This solution isn't a solution :

```
sage: [u[a].variables() for u in Solam]
[(), (a, c, d), (a, c, d), (a, c, d)]
```

The three last "solutions" are just rewritings of `eq3`

; the last is a triviality :

```
sage: Solam[0]
{a: 0}
```

leading to `{b:0} when substituting in the first equation, unacceptable as noted before...

Calling

```
sage: Solamp = solve(eq3, a, algorithm="maxima", solution_dict=True, to_poly_solve=True)
```

didnt't return after > 30 minutes

Similarly, `algorithm="sympy"`

didn't return after > 30 minutes. Trying `giac`

fails (Giac's solver seems limited to polynomials...). Frica's solver gives an answer seeminlgly using intermediary computations, which do not seem to be capturable "simply" by Sage (to be explored further...).

Another possibility is to use Sage's apparatus for polynomials :

```
sage: R1.<c0, c1>=QQbar[]
sage: R2.<c2, c3>=FractionField(R1)[]
sage: J=R2.ideal([R2(u.subs(dict(zip((a, b, c, d), (c0,c1, c2, c3))))) for u in (eq1, eq2)])
```

This ideal, which is, in some sense, the "best" solution to your problem, *can* be defined in Sage

```
sage: J
Ideal (c2 + (-c0^3 + (-3)*c0^2*c1 + (-4)*c0*c1^2)/(c0^2 + 3*c0*c1 + 2*c1^2), c3 + (-c0^2*c1 + (-2)*c0*c1^2 + (-2)*c1^3)/(c0^2 + 3*c0*c1 + 2*c1^2)) of Multivariate Polynomial Ring in c2, c3 over Fraction Field of Multivariate Polynomial Ring in c0, c1 over Algebraic Field
```

One can even be sure that this solution is a finite set :

```
sage: J.dimension()
verbose 0 (3998: multi_polynomial_ideal.py, groebner_basis) Warning: falling back to very slow toy implementation.
verbose 0 (1082: multi_polynomial_ideal.py, dimension) Warning: falling back to very slow toy implementation.
0
```

Alas :

```
sage: J.variety()
```

[ Snip... ]

```
NotImplementedError: root finding for this polynomial not implemented
```

OTOH, I'm far from an expert in this domain, and might have gone in completely wrong direction...

My tentative concusion isn't that Sage cannot solve this system, but rather that **I** can't solve this system with Sage but by using its `mathematica`

interface to the Wolfram engine...