kill(all)$
numer:true$ decimales:6$ prec:10^decimales$
a:1$ b:2$ N:10$
h:(b-a)/N$ 
t[0]:1$ y[0]:1$
define(f(t,y),0.7*y - t² + 1)$
y(t):=1.429*t^2+4.082*t-4.426*exp(.7*t)+4.402$
heun[0]:1$ E[0]:1$
print(" i ","y(t) Exacta", "  y(t) Euler","    y(t) Heun")$
for k:0 thru 10 do (
      t[k]:t[0]+k*h,
      t[k+1]:t[k]+h,
      y[k+1]:y[k]+h*f(t[k],y[k]),
      E[k+1]:heun[k]+h*f(t[k],heun[k]),
      heun[k+1]:heun[k]+(h/2)*(f(t[k],heun[k])+f(t[k+1],E[k+1])),
      print(t[k],float(round(y(t[k])*prec)/prec),"   ",
            float(round(y[k]*prec)/prec),"      ",
            float(round(heun[k]*prec)/prec))
)$
sol:makelist([t[j],heun[j]],j,0,N)$
print("Solución exacta versus solución Heun:")$
wxplot2d([[discrete,sol],y(t)],[t,a,b],[style,[points,2,2],
[lines,1,1]],[legend,"y(t) Heun","y(t) exacta"],
[xlabel,"t"], [ylabel,"y"]);