Taylor Method of order N

N:2$
f:(x+y-1)^2$
f:[f]$
for i:1 thru N-1 do
(
   g:expand(diff(f[i],x)+diff(f[i],y)*f[1]),
   f: append(f,[g])
)$
T: y$
kill(h)$
for i:1 thru N do
(
   T: T+h^i/i!*f[i]
)$
T;

Implement the Taylor method

x0: 0.0$
y0: 2.0$
h0: 0.05$
steps: 10$
for i:1 thru steps do
(
   y0: at( T, [x=x0, y=y0, h=h0] ),
   x0: x0+h0,
   print( "x = ", x0, "   y = ", y0 )
);