Equações diferenciais ordinárias: Método de Runge-Kutta de ordem 4

Código para solução de equações diferenciais ordinários com valor inicial usando o método de Runge-Kutta de quarta ordem.

import numpy as np

import math

def f(x,y):

    return x-2*y+1

x=0.0

y=1.0

h=0.1

n=10

p=np.zeros((n+1,3))

for i in range(1,n+1):

    k1=f(x,y)    

    k2=f(x+h/2,y+h*k1/2)

    k3=f(x+h/2,y+h*k2/2)

    k4=f(x+h,y+h*k3)    

    y=y+h*(k1+2*k2+2*k3+k4)/6

    x=h*i

    p[i,0]=i    

    p[i,1]=x

    p[i,2]=y

print(p)