Using the `solution_dict`

option of `solve`

makes it easier to use the solutions afterwards.

Here is an example based on the code in the question.

First we need to declare variables.

```
sage: var('B0 Q0 A0 P0 sigma0 B2 Q2 A2 P2 kf Q0d kp P0d C Q2d P2d')
```

Define equations:

```
sage: eqn1 = B0*Q0 == A0*P0 + C*(3*sigma0**2 + 1)
sage: eqn2 = B2*Q2 == A2*P2 + 2*C
sage: eqn3 = B0*kf*Q0d == A0*kp*P0d + 6*kp*C*sigma0
sage: eqn4 = B2*kf*Q2d == A2*kp*P2d
```

Solve and check out solutions:

```
sage: sol = solve([eqn1, eqn2, eqn3, eqn4], [B0, A0, B2, A2], solution_dict=True)
sage: len(sol)
1
sage: sol
[{B0: -(3*C*P0d*kp*sigma0^2 - 6*C*P0*kp*sigma0 + C*P0d*kp)/(P0*Q0d*kf - P0d*Q0*kp),
A0: -(3*C*Q0d*kf*sigma0^2 - 6*C*Q0*kp*sigma0 + C*Q0d*kf)/(P0*Q0d*kf - P0d*Q0*kp),
B2: -2*C*P2d*kp/(P2*Q2d*kf - P2d*Q2*kp),
A2: -2*C*Q2d*kf/(P2*Q2d*kf - P2d*Q2*kp)}]
```

Take the first solution (a dictionary):

```
sage: sola = sol[0]
sage: sola
{B0: -(3*C*P0d*kp*sigma0^2 - 6*C*P0*kp*sigma0 + C*P0d*kp)/(P0*Q0d*kf - P0d*Q0*kp),
A0: -(3*C*Q0d*kf*sigma0^2 - 6*C*Q0*kp*sigma0 + C*Q0d*kf)/(P0*Q0d*kf - P0d*Q0*kp),
B2: -2*C*P2d*kp/(P2*Q2d*kf - P2d*Q2*kp),
A2: -2*C*Q2d*kf/(P2*Q2d*kf - P2d*Q2*kp)}
```

Values of the solutions:

```
sage: sola[A0]
-(3*C*Q0d*kf*sigma0^2 - 6*C*Q0*kp*sigma0 + C*Q0d*kf)/(P0*Q0d*kf - P0d*Q0*kp)
sage: sola[B0]
-(3*C*P0d*kp*sigma0^2 - 6*C*P0*kp*sigma0 + C*P0d*kp)/(P0*Q0d*kf - P0d*Q0*kp)
sage: sola[A2]
-2*C*Q2d*kf/(P2*Q2d*kf - P2d*Q2*kp)
sage: sola[B2]
-2*C*P2d*kp/(P2*Q2d*kf - P2d*Q2*kp)
```

Simple substitutions:

```
sage: B = sola[B0].subs({P0d: 0})
sage: B
6*C*kp*sigma0/(Q0d*kf)
sage: A = sola[A0].subs({P0d: 0, P0: 1})
sage: A
-(3*C*Q0d*kf*sigma0^2 - 6*C*Q0*kp*sigma0 + C*Q0d*kf)/(Q0d*kf)
```

Noticing that one of the terms in `A`

looks like `B`

, we are tempted to substitute.

Unfortunately:

```
sage: A.subs({B: B0})
-(3*C*Q0d*kf*sigma0^2 - 6*C*Q0*kp*sigma0 + C*Q0d*kf)/(Q0d*kf)
```

Expanding and matching an exact term works:

```
sage: A.expand()
-3*C*sigma0^2 + 6*C*Q0*kp*sigma0/(Q0d*kf) - C
sage: A.expand().subs({B*Q0: B0*Q0})
-3*C*sigma0^2 + B0*Q0 - C
```