Tu script no me funcionaba, numacept hacia que accederia a ydibujo fuera de su rango.
Así que le agregue una condición para que parara una vez que excediera ese limite.
No le estabas mandado nada a plot, así que obviamente no ponía nada.
Imprimí los valores de h, x0 y xf y bueno h tiende a cero muy rápido lo que en 10000 iteraciones no iba a llevar a x0 a superar a xf.
Código Python:
Ver original# coding: utf-8
from pylab import *
# Insertamos la funcion con las condiciones oportunas
g = 9.81
l1 = 3
l2 = 7
m1 = 10000
m2 = 200
r = 0.75
c = 2 # medidas de m1
b = 1
mb = 50 # masa de la barra
def fcn(x, y):
return array ([y[1], -g*(m1*l1-m2*l2 + mb*(l1+l1)/2)*cos(y[0])/(m1*l1**2+m2*l2**2+(2/5)*m2*r**2+m1/12*(b**2+c**2) + 1/12*mb*(l1+l2)**2)])
# Definimos los coeficientes de la matriz RK
# c es un vector de dimension 7
c2 = 1.0E0 / 5.0E0
c3 = 3.0E0 / 10.0E0
c4 = 4.0E0 / 5.0E0
c5 = 8.0E0 / 9.0E0
c6 = 1.0E0
c7 = 1.0E0
# La matriz a es de dimension 7
a21 = 1.0E0 / 5.0E0
a31 = 3.0E0 / 40.0E0
a41 = 44.0E0 / 45.0E0
a51 = 19372.0E0 / 6561.0E0
a61 = 9017.0E0 / 3168.0E0
a71 = 35.0E0 / 384.0E0
a32 = 9.0E0 / 40.0E0
a42 = -56.0E0 / 15.0E0
a52 = -25360.0E0 / 2187.0E0
a62 = -355.0E0 / 33.0E0
a43 = 32.0E0/9.0E0
a53 = 64448.0E0/6561.0E0
a63 = 46732.0E0/5247.0E0
a73 = 500.0E0/1113.0E0
a54 = -212.0E0/729.0E0
a64 = 49.0E0/176.0E0
a74 = 125.0E0/192.0E0
a65 = -5103.0E0/18656.0E0
a75 = -2187.0E0/6784.0E0
a76 = 11.0E0/84.0E0
# el vector b es el del paso mayor (orden5)
b1 = 35.0E0 / 384.0E0
b3 = 500.0E0 / 1113.0E0
b4 = 125.0E0 / 192.0E0
b5 = -2187.0E0 / 6784.0E0
b6 = 11.0E0 / 84.0E0
# El vector be es el de orden menor (orden 4)
be1 = 5179.0E0 / 57600.0E0
be2 = 0
be3 = 7571.0E0 / 16695.0E0
be4 = 393.0E0 / 640.0E0
be5 = -92097.0E0 / 339200.0E0
be6 = 187.0E0 / 2100.0E0
be7 = 1.0E0 / 40.0E0
numrecha = 0
numacept = 0
# intervalo en el que hay que realizar los calculos
x0 = 0.0
xf = 5.0
# errores abs y rel
EABS = 1e-8
EREL = 1e-9
y0 = array([arccos(-1.0)/4, 0.0])
y = zeros([size(y0)])
n = 1000
tol = EABS + EREL * norm(y0)
g1 = fcn(x0, y0)
h = min(0.5, (tol/(1e-38+norm(g1)))**(1/5))
print h
xab = zeros([n, 1])
yab = zeros([n, 2])
zab = zeros([n, 2])
fallo = False
ydibujo = zeros([10000, 2])
while x0 < xf and numacept < 10000:
g2 = fcn(x0+c2*h, y0+h*a21*g1)
g3 = fcn(x0+c3*h, y0+h*(a31*g1+a32*g2))
g4 = fcn(x0+c4*h, y0+h*(a41*g1+a42*g2+a43*g3))
g5 = fcn(x0+c5*h, y0+h*(a51*g1+a52*g2+a53*g3+a54*g4))
g6 = fcn(x0+c6*h, y0+h*(a61*g1+a62*g2+a63*g3+a64*g4+a65*g5))
g7 = fcn(x0+c7*h, y0+h*(a71*g1+a73*g3+a74*g4+a75*g5+a76*g6))
# Calculamos la y
y5 = y0+h*(b1*g1+b3*g3+b4*g4+b5*g5+b6*g6)
y4 = y0+h*(be1*g1+be2*g2+be3*g3+be4*g4+be5*g5+be6*g6+be7*g7)
# tolerancia y estimación
est = norm(y5-y4)
tol = EABS+EREL*max(1, norm(y5))
if est <= tol:
# AQUI EL PASO SE ACEPTA
x0 = x0 + h
y0 = y5
ydibujo[numacept, :] = y5
g1 = g7
fac = min(2, 0.9*(tol/(est+1e-38))**(1/5))
if fallo:
fac = min(1, fac)
h = fac*h
if x0+h > xf:
h = xf-x0
fallo = False
numacept = numacept + 1
else:
# El paso ha sido rechazado por tanto hay que recalcular el paso.
h = max(0.1, 0.9*(tol/est)**(1/5))*h
if fallo:
print "Dos pasos rechazados consecutivos"
fallo = True
numrecha = numrecha+1
print numrecha
print numacept
plot(ydibujo)
show()
También encontré una implementación usando funciones lambda.
https://rosettacode.org/wiki/Runge-Kutta_method#Python