1001.sagews

b1, b0, N, B, p = var('b1 b0 N B p')
eqB = B == b1 - ((b1 - b0) / (p*N + 1-p))
print eqB

B == b1 + (b0 - b1)/(N*p - p + 1)

Q, q0, q1 = var('Q q0 q1')
sol = solve([eqB, B==Q, b0==q0, b1==q1, q0==10, q1==40, p==0.35], Q, q0, q1, B, b0, b1, p)[0]
eqQ = sol[0]
eqq0 = sol[1]
print eqQ
print eqq0
plot(eqQ.rhs(), (N, 1, 50))

Q == 40*(7*N - 2)/(7*N + 13) q0 == 10

i2, i0, i1 = var('i2 i0 i1')
sol = solve([eqB, B==i2, b0==i0, b1==i1, i0==8, i1==3, p==0.35], i2, i0, i1, B, b0, b1, p)[0]
eqi2 = sol[0]
eqi0 = sol[1]
print eqi2
print eqi0
plot(eqi2.rhs(), (N, 1, 50))

i2 == (21*N + 139)/(7*N + 13) i0 == 8

c, ce, RE = var('c ce RE')
eqRE = RE == (q0*(1-c) + Q*c)/(i0*(1-c) + i2*c + c*(i2/ce)*(N*(N+1)/2))
print eqRE

RE == -2*((c - 1)*q0 - Q*c)/((N + 1)*N*c*i2/ce - 2*(c - 1)*i0 + 2*c*i2)

RE0 = var('RE0')
sol = solve([eqRE, RE0==RE, eqQ, eqq0, eqi2, eqi0, c==0.5, ce==40], RE0, RE, Q, q0, i2, i0, c, ce)
eqRE0 = sol[0][0]
print eqRE0
plot(eqRE0.rhs(), (N, 1, 30))

RE0 == 4000*(7*N + 1)/(21*N^3 + 160*N^2 + 6299*N + 19440)