function a = f(x, y)
    a = 3 * x - y + 8
endfunction
 
function a = k1(x, y, h)
    a = h * f(x, y)
endfunction
 
function a = k2(x, y, h)
    a = h * f(x + h / 2, y + k1(x, y, h) / 2)
endfunction
 
function a = k3(x, y, h)
    a = h * f(x + h / 2, y + k2(x, y, h) / 2)
endfunction
 
function a = k4(x, y, h)
    a = h * f(x + h, y + k3(x, y, h))
endfunction
 
 
x = 0:0.1:1
y(1) = 3
h = 0.1
 
for i = 1:10
    y(i + 1) = y(i) + (k1(x(i), y(i)) + 2 * k2(x(i), y(i)) + 2 * k3(x(i), y(i)) + k4(x(i), y(i))) / 6
    disp([x(i), y(i)])
end