In the past, the price and quantity are both equal in the demand and supply equations, which makes our computation fairly easy. In this case, we actually have two different prices and the same quantity. This difference requires us to create a new function that solves a system of equations. Let’s walk step by step through the creation of this function.

import sympy
p = sympy.Symbol("p")
q = sympy.Symbol("q")

def eqSolve(eq1,eq2,tax):
    demandP = sympy.solve(eq1-q,p)[0]
    supplyP = sympy.solve(eq2-q,p)[0]
    print(demandP)
    print(supplyP)
eqSolve(10-p,2*p,2)

Our first equation is going to be our demand equation, the second will be our supply equation. The first two lines are taking the left hand side variable, quantity, and moving it to the right side. The equation for demand we are using, for example, is q = 10-p but for sympy we want it in the form 10-p-q when we solve.

The second argument, p, is telling sympy that we want to solve this equation for the variable p. While our equation tells us what q equals, we are going to need what the price equals for both demand and supply.

Now let’s add on the tax. We are going to specify the tax in this way. The equilibrium price without tax is p. The price the consumer pays is p+tax

_consumer

, and the price the producer receives is p-tax

_producer

. The tax is added for the consumer because they pay more, where as for the producer they get less revenue so it is subtracted. Finally, the two taxes must add up to the total tax so we get: tax = tax

_consumer

+ tax

_producer

.

Right now we have p = -q + 10 and p = q/2, this becomes p+tax

_consumer

= -q + 10 and p + tax

_producer

= q/2. We need to move the taxes to the right in our code, as shown below.

import sympy
p = sympy.Symbol("p")
q = sympy.Symbol("q")
cTax = sympy.Symbol("cTax")
pTax = sympy.Symbol("pTax")

def eqSolve(eq1,eq2,tax):
    demandP = sympy.solve(eq1-q,p)[0]
    supplyP = sympy.solve(eq2-q,p)[0]
    demandP = demandP-cTax
    supplyP = supplyP+pTax
    print(demandP)
    print(supplyP)
eqSolve(10-p,2*p,2)