clear;

x0 = input("initial value of x = "); 
y0 = input("initial value of y = "); 
x1 = input("upper limit of x = ");
N  = input("number of steps = ");


function z = f1(x, y)
    z = x*y;
end

x = zeros(N, 1);
y = zeros(N, 1);

h = (x1-x0) / N;
x(1) = x0;
y(1) = y0;

for i=1:N
    k1 = h * f1(x(i),y(i));  
    k2 = h * f1(x(i) + h/2.0, y(i) + 0.5 * k1);  
    k3 = h * f1(x(i) + h/2.0, y(i) + 0.5 * k2);  
    k4 = h * f1(x(i) + h, y(i) + k3);  
    y(i+1) = y(i) + (k1 + 2.0*k2 +2.0*k3 + k4 ) / 6.0;
    x(i+1) = x(i) + h;
end

[x y]