f(x,y):=x+y;
x[0]:0;y[0]:1;

a(x,y):=''(diff(f(x,y),x,1));

b(x,y):=''(diff(f(x,y),y,1));

c(x,y):=''(diff(f(x,y),x,2));

d(x,y):=''(diff(f(x,y),y,2));

e(x,y):=''(diff(a(x,y),y,1));

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

yprime(x[0], y[0]);

yyprime(x,y):=a(x,y)+b(x,y)*f(x,y);

yyprime(x[0], y[0]);

ytaylor(x):= y[0] + yprime[0]*(x-x[0]) + yyprime(x[0], y[0])*(x-x[0])^2 / 2!;

ytaylor(0.5);

h: float(1/3);

for n: 1 thru 3 do (
x[n]: x[n-1] + h,
y[n]: y[n-1] + h + yprime(x[n-1], y[n-1])*h + yyprime(x[n-1], y[n-1])*h^2 / 2!
);

makelist([x[n], y[n]], n, 0, 3);

ode2('diff(y, x) = f(x,y), y, x);

ic1(%, x = x[0], y = y[0]);