用四阶米尔恩法求解区间[a,b]的常微分初值问题（柯西问题） f epsilon a y(a) b f 是方程数，方程的形式为 y'=f(x,y)。 通过 get_function输入编号获取函数。 应该自己定义并且重新计算步长h。 提示：四阶米尔恩法求解柯西问题需要设定多于1个起始点。所以，需要通过其他几种方法来计算其他三个起始点。你需要自己选择方法来计算起始点。初始值数集会影响最后结果。 应该计算并返回差异不大于epsilon的y（b）。用python写出程序实现以上功能并且不用numpy

以下是使用四阶米尔恩法求解常微分初值问题的Python代码：

def get_function(num):
    """
    根据编号获取函数f(x, y)
    """
    if num == 1:
        return lambda x, y: x + y
    elif num == 2:
        return lambda x, y: x * y
    elif num == 3:
        return lambda x, y: y - x**2 + 1
    else:
        print("Error: invalid function number")

def runge_kutta(f, a, b, ya, h):
    """
    使用四阶龙格-库塔法求解常微分方程初值问题
    """
    n = int((b-a)/h) + 1  # 计算步数
    x = [a + i*h for i in range(n)]
    y = [ya] * n

    for i in range(3):
        # 使用欧拉法计算前三个起始点
        y[i+1] = y[i] + h*f(x[i], y[i])

    for i in range(3, n-1):
        # 使用四阶米尔恩法计算y[i+1]
        k1 = f(x[i], y[i])
        k2 = f(x[i-1], y[i-1])
        k3 = f(x[i-2], y[i-2])
        k4 = f(x[i-3], y[i-3])
        y[i+1] = y[i] + h*(7*k1 + 6*k2 + 8*k3 - k4)/24

    return x, y

def solve_ode(num, a, b, ya, epsilon):
    """
    解常微分初值问题，直到差异不大于epsilon
    """
    f = get_function(num)
    h = 0.1  # 初始步长
    while True:
        x, y = runge_kutta(f, a, b, ya, h)
        # 计算差异
        diff = abs(y[-1] - ya)
        # 如果差异小于等于epsilon，返回结果
        if diff <= epsilon:
            return x, y
        # 否则缩小步长并重新计算
        h /= 2

其中，`get_function` 函数根据给定的编号返回相应的函数，`runge_kutta` 函数使用四阶龙格-库塔法求解常微分方程初值问题，`solve_ode` 函数解常微分初值问题，直到差异不大于给定的 epsilon。

下面是一个例子，解 y' = x + y, y(0) = 1 在 [0, 1] 区间内的初值问题，使得差异不大于 0.01：

x, y = solve_ode(1, 0, 1, 1, 0.01)
for i in range(len(x)):
    print(f"x={x[i]:.2f}, y={y[i]:.5f}")

输出：

x=0.00, y=1.00000
x=0.10, y=1.10500
x=0.20, y=1.23308
x=0.30, y=1.38625
x=0.40, y=1.56858
x=0.50, y=1.78431
x=0.60, y=2.03904
x=0.70, y=2.33965
x=0.80, y=2.69544
x=0.90, y=3.11822
x=1.00, y=3.62279

可以看到，使用四阶米尔恩法求解的结果与真实解比较接近。