Code:

p_string = 'FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFE FFFFFC2F'
p = ZZ( '0x' + p_string.replace(' ', ''))
E = EllipticCurve( GF(p), [0, 7] )
n_string = 'FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFE BAAEDCE6 AF48A03B BFD25E8C D0364141'
n = ZZ( '0x' + n_string.replace(' ', '') )

xP = 0x79be667ef9dcbbac55a06295ce870b07029bfcdb2dce28d959f2815b16f81798
yP = 0x483ada7726a3c4655da4fbfc0e1108a8fd17b448a68554199c47d08ffb10d4b8
G = E.point( (xP, yP) )

def egcd(a, b):
"Euclidean greatest common divisor"
if a == 0:
return (b, 0, 1)
else:
g, y, x = egcd(b % a, a)
return (g, x - (b // a) * y, y)

def modinv(a, m):
"Modular inverse"
# in Python 3.8 you can simply return pow(a,-1,m)
g, x, y = egcd(a, m)
if g != 1:
raise Exception('modular inverse does not exist')
else:
return x % m

#making point from two coordinates as x,y to P
px1 = 0x79be667ef9dcbbac55a06295ce870b07029bfcdb2dce28d959f2815b16f81798
py1 = 0x483ada7726a3c4655da4fbfc0e1108a8fd17b448a68554199c47d08ffb10d4b8

P = E.point( (px1,py1) )

px2 = 0xc6047f9441ed7d6d3045406e95c07cd85c778e4b8cef3ca7abac09b95c709ee5
py2 = 0x1ae168fea63dc339a3c58419466ceaeef7f632653266d0e1236431a950cfe52a

P2 = E.point( (px2,py2) )

#add Two point toget P3:

P3= P+P2

#substract
P3= P-P2

# div by scalar

P4= P2*modinv(2,n)
# 2*(1/2)=1
assert P4==P

#mul by scalar

P5=P1*2
assert P5==P2