import numpy as np
import pylab as pl

Y, X = np.mgrid[-10:10:100j, -10:10:100j]
# Domain Parameters
a = 3
bb = 7
def f(x):
    return x**3 -a*x + bb

###first point P
px = -2.0
py = np.sqrt(f(px))
###second point Q
qx = -0.5
qy = np.sqrt(f(qx))

print "point p x:%s y:%s" % (px,py)
print "point q x:%s y:%s" % (qx,qy)
### Slope
k = (qy - py)/(qx - px)
b = -px*k + py
print "Slope: %s" % k

##Point Addition
realx = pow(k,2) - qx - px
realy = qy + k * (realx-qx)

## Plotting
pl.contour(X, Y, Y**2 - f(X), levels=[0])
plt.plot(px,py, linestyle='-', marker='o', color='y')
plt.plot(qx,qy, linestyle='-', marker='o', color='r')
plt.plot([realx,realx],[-realy,realy], linestyle='--', marker='o', color='r')

pl.plot([px,realx],[py,realy],linestyle="--")
pl.show()