kill(all)$
numer:true$ decimales:6$ prec:10^decimales$
a:0$ b:1$ N:10$
h:(b-a)/N$ 
t[0]:0$ y[0]:1$
y(t):=.5*exp(-2*t)+.5*exp(-4*t)$
define(f0(t,y),exp(-2*t)-4*y);
define(f1(t,y),diff(f0(t,y),t,1)+ diff(f0(t,y),y,1)*f0(t,y));
define(f2(t,y),diff(f1(t,y),t,1)+ diff(f1(t,y),y,1)*f0(t,y));
define(f3(t,y),diff(f2(t,y),t,1)+ diff(f2(t,y),y,1)*f0(t,y));
print("i   ","   Analítica","    Taylor")$
for k:0 thru N do (
      t[k]:t[0]+k*h,
      t[k+1]:t[0]+h,
      y[k+1]:y[k]+h*(f0(t[k],y[k])
             +(h/2)*f1(t[k],y[k])
             +(h^2/6)*f2(t[k],y[k])
             +(h^3/24)*f3(t[k],y[k])),
      print(round(t[k]*10)/10,"    ",
            float(round(y(t[k])*prec)/prec),"   ",
            float(round(y[k]*prec)/prec))        
)$

wxdraw2d(
color=red,
key="Solución Taylor",
point_type=filled_circle,
point_size=1,
points(makelist([t[j],y[j]],j,0,N)),
color=blue,
key="Solución Exacta",
explicit(y(t),t,a,b),
grid=true,
xlabel="t",
ylabel="y",
title="Solución analítica versus solución numérica"
);