Second order Taylor method for solving ODE

Question

Can anyone help me modify this code to second-order Taylor method?

# Python Code to find the approximation
# of an ordinary differential equation
# using Euler method.

# Consider a differential equation
# dy / dx =(x + y + xy)
def func(x, y):
    return (x + y + x * y)


# Function for euler formula
def euler(x0, y, h, x):
    temp = -0

    # Iterating till the point at which we
    # need approximation
    while x0 < x:
        temp = y
        y = y + h * func(x0, y)
        x0 = x0 + h

    # Printing approximation
    print("Approximate solution at x = ", x, " is ", "%.6f" % y)


# Driver Code
# Initial Values
x0 = 0
y0 = 1
h = 0.025

# Value of x at which we need approximation
x = 0.1

euler(x0, y0, h, x)


temp=-0
def second_order(x0,y,h,x):
    while x0 < x:
        temp = y
        y = y + h * func(x0, y)+(h^2/2)*func(x0,y)
        x0 = x0 + h
    print("Approximate solution at x = ", x, " is ", "%.6f" % y)

second_order(x0,y0,h,x)

Answer 1

Use h**2 or pow(h,2) for $h^2$ as, ^ is an EXOR operator in Python and is meant for integers

After doing that and running the same code I got

Approximate solution at x =  0.1  is  1.111673
Approximate solution at x =  0.1  is  1.113123

Answer 2

Taylor approximation of order 2 looks like this $$f(x) \approx f(x_0) + f'(x_0)(x-x_0) + f''(x_0)(x-x_0)^2/2$$

At $x= x_0+h$ this reads $$f(x_0+h) \approx f(x_0) + hf'(x_0) + \frac{h^2}{2}f''(x_0)$$

Thus you need to calculate $f''$ (your formula use $f'$), which you can do using the differential equation.

That is, you would like to use

def func2(x, y):
    return 1+y + (1+x)* func(x,y)

(Also, the syntax $h^2$ should be replaced by python syntax $h*h$).