taylor.py

from sympy import *
from sympy.plotting import plot3d
import mpmath as mp
import numpy as np
x = Symbol('x')
y = Symbol('y')
f = exp(x / 2) * sin(y)

fx = f.diff(x)
fy = f.diff(y)
fxx = fx.diff(x)
fyy = fy.diff(y)
fxy = fx.diff(y)

print("f = ", f)
print("fx = ", fx)
print("fy = ", fy)
print("fxx = ", fxx)
print("fyy = ", fyy)
print("fxy = ", fxy)

x0 = 1
y0 = 1

# Callable function
g = Lambda((x, y), f)
gx = Lambda((x, y), fx)
gy = Lambda((x, y), fy)
gxx = Lambda((x, y), fxx)
gyy = Lambda((x, y), fyy)
gxy = Lambda((x, y), fxy)
def q(x, y):
    return g(x0, y0) + gx(x0, y0) * (x - x0) + gy(x0, y0) * (y - y0) + gxy(x0, y0) * (x - x0) * (y - y0) + 1/2 * gxx(x0, y0) * (x - x0)**2 + 1/2 * gyy(x0, y0) * (y - y0)**2

qs = q(x, y)
approx= plot3d(qs)