1216.sagews

def runge_kutta_four(a,b,N,alpha):
    #create functions we will need
    f = lambda x,y: 1.0 + (y/x) + (y/x)^2 #this is the first equation
    g = lambda x: x*tan(ln(x)) # this is the exact solution for the IVP
    #continuing with page 279
    h = (b-a)/N
    t = a
    omega = alpha
    #make a list to store values
    w = []
    w.append((0,omega, n(abs(g(1)-omega)))) # this stores y(1) = 0 & actual error
    for i in range(1,N):
        K1 = h * f(t,omega)
        K2 = h * f(t + (h/2), omega + (K1/2))
        K3 = h * f(t + (h/2), omega + (K2/2))
        K4 = h * f(t + h, omega + K3)
        omega = omega + (K1 + K2 + K3 + K4)/6 # compute omega_i
        w.append(( i, omega, n(abs(g(1)-omega)) ) ) #store values
        t = a + i * h # compute t_i
    return w

ans = runge_kutta_four(1,4,300,0)
mat_ans = matrix(ans)
latex_mat = latex(mat_ans)