direction field

h = 5 (slider from 0.1 to 7, step 0.1)
a = 4
tau = 5
f(x, y) = 1 / tau * (-y + a * ℯ^((-x) / tau))
y0 = 1
x_0: (0, y0)
x_1: x_0 + (h, h * f(x_0))
x_2: x_1 + h * (1, f(x_1))
x_3: x_2 + h * (1, f(x_2))
x_4: x_3 + h * (1, f(x_3))
x_5: x_4 + h * (1, f(x_4))
p6: x_5 + h * (1, f(x_5))
p7: p6 + h * (1, f(p6))
p8: p7 + h * (1, f(p7))
p9: p8 + h * (1, f(p8))
x_{10}: p9 + h * (1, f(p9))
p11: x_{10} + h * (1, f(x_{10}))
p12: p11 + h * (1, f(p11))
p13: p12 + h * (1, f(p12))
p14: p13 + h * (1, f(p13))
x_{15}: p14 + h * (1, f(p14))
p16: x_{15} + h * (1, f(x_{15}))
p17: p16 + h * (1, f(p16))
p18: p17 + h * (1, f(p17))
p19: p18 + h * (1, f(p18))
x_{20}: p19 + h * (1, f(p19))
p21: x_{20} + h * (1, f(x_{20}))
p22: p21 + h * (1, f(p21))
p23: p22 + h * (1, f(p22))
p24: p23 + h * (1, f(p23))
x_{25}: p24 + h * (1, f(p24))
p26: x_{25} + h * (1, f(x_{25}))
p27: p26 + h * (1, f(p26))
p28: p27 + h * (1, f(p27))
p29: p28 + h * (1, f(p28))
x_{30}: p29 + h * (1, f(p29))
p31: x_{30} + h * (1, f(x_{30}))
p32: p31 + h * (1, f(p31))
p33: p32 + h * (1, f(p32))
p34: p33 + h * (1, f(p33))
x_{35}: p34 + h * (1, f(p34))
p36: x_{35} + h * (1, f(x_{35}))
exact = Function((y0 + a / tau * x) * ℯ^((-x) / tau), 0, Infinity)
dir0 = Vector(x_0, x_0 + (h, h * f(x_0)))
dir1 = Vector(x_1, x_1 + h * (1, f(x_1)))
u = Vector(x_2, x_3)
v = Vector(x_3, x_4)
w = Vector(x_4, x_5)
z = Vector(x_5, p6)
b = Vector(p6, p7)
c = Vector(p7, p8)
d = Vector(p8, p9)
e = Vector(p9, x_{10})
g = Vector(x_{10}, p11)
i = Vector(p11, p12)
j = Vector(p12, p13)
k = Vector(p13, p14)
l = Vector(p14, x_{15})
m = Vector(x_{15}, p16)
n = Vector(p16, p17)
p = Vector(p17, p18)
q = Vector(p18, p19)
r = Vector(p19, x_{20})
s = Vector(x_{20}, p21)
t = Vector(p21, p22)
u_1 = Vector(p22, p23)
v_1 = Vector(p23, p24)
w_1 = Vector(p24, x_{25})
z_1 = Vector(x_{25}, p26)
a_1 = Vector(p26, p27)
b_1 = Vector(p27, p28)
c_1 = Vector(p28, p29)
d_1 = Vector(p29, x_{30})
e_1 = Vector(x_{30}, p31)
f_1 = Vector(p31, p32)
p37: p36 + h * (1, f(p36))
p38: p37 + h * (1, f(p37))
p39: p38 + h * (1, f(p38))
x_{40}: p39 + h * (1, f(p39))
p41: x_{40} + h * (1, f(x_{40}))
p42: p41 + h * (1, f(p41))
p43: p42 + h * (1, f(p42))
p44: p43 + h * (1, f(p43))
x_{45}: p44 + h * (1, f(p44))
p46: x_{45} + h * (1, f(x_{45}))
p47: p46 + h * (1, f(p46))
p48: p47 + h * (1, f(p47))
g_1 = Vector(p32, p33)
h_1 = Vector(p33, p34)
i_1 = Vector(p34, x_{35})
k_1 = Vector(p37, p38)
l_1 = Vector(p38, p39)
m_1 = Vector(p39, x_{40})
n_1 = Vector(x_{40}, p41)
p_1 = Vector(p41, p42)
q_1 = Vector(p42, p43)
r_1 = Vector(p43, p44)
s_1 = Vector(p44, x_{45})
t_1 = Vector(x_{45}, p46)
u_2 = Vector(p46, p47)
v_2 = Vector(p47, p48)
p49: p48 + h * (1, f(p48))
x_{50}: p49 + h * (1, f(p49))
w_2 = Vector(p48, p49)
z_2 = Vector(p49, x_{50})
j_1 = Vector(x_{35}, p36)
a_2 = Vector(p36, p37)
dx = 1
dy = 0.25
o = false
dirfield = Sequence(Sequence[Segment[(a, b), (a + dx / 3, b + dx / 3 * f(a, b))], a, -4, 42, dx], b, -1.75, 4.5, dy)