f(x,y):= x+y;

x[0]: 0; y[0]: 1; h: 0.25;

for n: 1 thru 4 do (
R1: h*f(x[n-1], y[n-1]),
R2: h*f(x[n-1] + .5*h, y[n-1]+.5*R1),
R3: h*f(x[n-1] + .5*h, y[n-1]+.5*R2),
R4: h*f(x[n-1] + h, y[n-1]+R3),
y[n]: y[n-1] + 1/6 * (R1 + 2*R2 + 2*R3 + R4),
x[n]: x[n-1] + h
);

data: makelist([x[n], y[n]], n, 0, 4);