Method to find solution of ordinary differential equations

ODE-Heun.py

#Heun Method for solutions approach
import math, sys 
import numpy as np

def Heun(f, y0, T, n):
    """Solve y'=f(t,y), y(0)=y0, with n steps until t=T."""
    t = np.zeros(n+1)
    y = np.zeros(n+1)  # y[k] is the solution at time t[k]
    dt = T/float(n)
    y[0] = y0
    t[0] = 0
    for k in range(n):
        t[k+1] = t[k] + dt
        y[k+1] = y[k] + (dt/2)*(f(y[k],t[k]) + f(y[k],t[k+1] + dt*f(y[k],t[k])))
    return y, t