1
var('xx yy x y')
2
   
3
@interact
4
def _(F=cos(x)*cos(y),x0=(0,1.5), y0=(0, 1.5),
5
      order=[1..10]):
6

7
    plotF = plot3d(F, (-2, 2), (-2, 2), adaptive=True, color='blue')
8
    F0 = float(F.subs(x=x0).subs(y=y0))
9
    P = (x0, y0, F0)
10
    dot = point3d(P, size=15, color='yellow')
11
    plot = dot + plotF
12
    approx = F0
13
    for n in range(1, order+1):
14
        for i in range(n+1):
15
            if i == 0:
16
                deriv = F.diff(y, n)
17
            elif i == n:
18
                deriv = F.diff(x, n)
19
            else:
20
                deriv = F.diff(x, i).diff(y, n-i)
21
            deriv = float(deriv.subs(x=x0).subs(y=y0))
22
            coeff = binomial(n, i)/factorial(n)
23
            approx += coeff * deriv * (x-x0)^i * (y-y0)^(n-i)
24
    plot += plot3d(approx, (x, -2, 2),
25
             (y, -2, 2), color='red', opacity=0.7)
26
    pretty_print("F(x,y) =",F)
27
    show(plot)
28

29