def noisyrk4(s,t,tau,derivsRK,i,vN0,vN1,vN2):
    """modified RK4 integrator including noise
    DEPENDENCIES
        none
    INPUTS
        s - initial state vector [delta]
        t - time
        tau - time step size
        derivsRK - RHS of ODE, fn defined somewhere
        i - bias current
        vN0,vN1,vN2 - voltage noise values at t=t, t=t+tau/2, t=t+tau
    OUTPUTS
        sout - new state vector [delta]"""
    
    half_tau = 0.5*tau
    
    F1 = derivsRK(s,t,i,vN0)                    # use current voltage noise
    t_half = t + half_tau
    stemp = s + half_tau*F1
    
    F2 = derivsRK(stemp,t_half,i,vN1)           # use half-tau step voltage noise
    stemp = s + half_tau*F2
    
    F3 = derivsRK(stemp,t_half,i,vN1)           # use half-tau step voltage noise
    t_full = t + tau
    stemp = s + tau*F3
    
    F4 = derivsRK(stemp,t_full,i,vN2)           # use full-tau step voltage noise
    sout = s + tau/6.*(F1 + F4 + 2.*(F2 + F3))
    return sout