# Verfahren von Heun zur numerischen Lösung von Anfangswertaufgaben.
#
#def heun(xk,yk,h,ys):
#    xk1 = xk + h
#    yk1p = yk + h * ys(xk,yk)
#    yk1 = 0.5*(yk+yk1p+h*ys(xk1,yk1p))
#    return (xk1,yk1p)
#
#def run_heun(K,h,x0,y0,ys):
#    if not isinstance(ys,list): ys = [ys]
#    if not isinstance(y0,list): y0 = [y0]
#    xk = x0
#    yk = y0
#    xks = x0
#    yks = [ [v] for v in y0 ]
#
#    for i in range(K):
#        for j in range(len(ys)):
#            xk,yk[j] = heun(xk,yk[j],h,ys[i])
#            xks.append(xk)
#            yks[j].append(yk[j])
#

def run_heun(K,h,x0,y0,ys):
    if not isinstance(y0,list): y0 = [y0]
    if not isinstance(ys,list): ys = [ys]

    xk = x0
    yk = y0
    xka = 0
    yka = [0]*len(yk)
    ykp = [0]*len(yk)

    xks = [x0]
    yks = [ [v] for v in y0 ]

    for i in range(K):
        xka = xk + h
        xks.append(xka)
        for j in range(len(ys)):
            yka[j] = yk[j] + h * ys[j](xk,*yk)
        for j in range(len(ys)):
            ykp[j] = 0.5*(yk[j] + yka[j] + h * ys[j](xka,*yka))
            yks[j].append(ykp[j])
        yk = ykp
        xk = xka

    return [xks] + yks