kill(all)$
numer:true$
a:1$ b:2$ N:10$
h:(b-a)/N$ 
t[0]:1$ 
define(f(t,y),0.7*y - t² + 1)$
y(t):=1.43*t^2+4.08*t-4.425*exp(.7*t)+4.4$
y[0]:1$
for k:1 thru 10 do (
      t[k]:t[0]+k*h,
     y[k]:y[k-1]+h*f(t[k-1],
     y[k-1])
)$
print("Salida numérica: ",
sol:makelist([t[j],y[j]],j,0,N))$
print("Solución exacta versus solución aproximada:")$
wxplot2d([[discrete,sol],y(t)],[t,a,b],[style,[points,2,2],
[lines,1,1]],[legend,"y(t) aproximada","y(t) exacta"],
[xlabel,"t"], [ylabel,"y"]);