It seems like this is just ridiculously slow, but I'm not sure why. I tried tweaking your code in a variety of ways and nothing helped. To try to figure out what was going on I wrote the following Runge-Kutta 4th order code in Cython; its pretty hackish but it makes it clear that the solution isn't doing anything very strange. The version below is hard-coded for your ll=2.

```
%cython
cdef double ll = 2.0
cpdef double leg_f1(double t, double y, double yp):
return yp
cpdef double leg_f2(double t, double y, double yp):
return ((ll**2 + ll)*y - 2.0*yp*t)/(t**2 - 1.0)
cpdef RK4_2d(f1, f2, double t_start, double y_start, double yp_start, double t_end, int steps):
cdef double step_size = (t_end - t_start)/steps
cdef double t_current = t_start
cdef double y_current = y_start
cdef double yp_current = yp_start
cdef list answer_table = []
cdef int j
answer_table.append([t_current,y_current, yp_current])
for j in range(0,steps):
k1=f1(t_current, y_current, yp_current)
k2=f1(t_current+step_size/2, y_current + k1*step_size/2, yp_current + k1*step_size/2)
k3=f1(t_current+step_size/2, y_current + k2*step_size/2, yp_current + k1*step_size/2)
k4=f1(t_current+step_size, y_current + k3*step_size, yp_current + k1*step_size/2)
y_current = y_current + (step_size/6)*(k1+2*k2+2*k3+k4)
k1=f2(t_current, y_current, yp_current)
k2=f2(t_current+step_size/2, y_current + k1*step_size/2, yp_current + k1*step_size/2)
k3=f2(t_current+step_size/2, y_current + k2*step_size/2, yp_current + k1*step_size/2)
k4=f2(t_current+step_size, y_current + k3*step_size, yp_current + k1*step_size/2)
yp_current = yp_current + (step_size/6)*(k1+2*k2+2*k3+k4)
t_current += step_size
answer_table.append([t_current,y_current,yp_current])
return answer_table
```

and then you can get a plot with:

```
ans = RK4_2d(leg_f1, leg_f2, 0.0, 0.0, 1.0, 1.0, 40)
list_plot([[x[0],x[1]] for x in ans], plotjoined = True)
```

Anyway for me it seems easier to implement such things directly than to wrestle with the GSL wrapping.

Just to check: Are you solving Legendre's DE?

yea, im just trying to learn Sage solving problems like this. for Pl(x)=x i found the eigenvalue l= 0.994838832387 using desolve_system_rk4. I just wanna to solve the problem using ode_solver and odeint.

how i do backtrace? the script doesn't crash.... it just calculate forever

Well... it seems to be that when the solver evaluates the function (eqa(t=1)) the solution diverges because the 1/t factor. But this doesn't happens with the desolver_system_rk4. i dont like this, maybe ode_solver is faster but.. it really has to have this kind of problems?