3739.sagews

import numpy as np

def f(t,y):
 return (y)

def RK4 (a,b,N,ya):
    y= np.empty(N+1)
    h=(b-a)/N
    y[0]=ya
    for i in range(N):
        t=a+i*h
        tdemi=t+h/2
        r1=f(t,y[i])
        r2=f(tdemi,y[i]+(h/2)*r1) 
        r3=f(tdemi,y[i]+(h/2)*r2)...

        r4=f(t+h,y[i]+h*r3)
        y[i+1]=y[i]+h/2*(r1+2*r2+2*r3+r4)
    return (y)

Y=RK4(0,1,10,1)

list_plot([(0.1*1,Y[1]) for 1 in range (11)])...