/* y'=x*y */
rhs:x*y;

/* exact solution (separation of variables) and check */
es:K*exp(x^2/2);
diff(es,x)-subst(y=es,rhs);

/* number of steps and increment dx */
n:80;
dx:0.05;

/* initialize arrays: y(0)=1 */
x[0]:0$
y[0]:1$
h[0]:y[0]$
rk[0]:y[0]$

for i:1 thru n do (
    
    x[i]:x[i-1]+dx,

    /* forward Euler */
    y[i]:y[i-1]+subst([x=x[i-1],y=y[i-1]],rhs)*dx,

    /* Heun */
    m1:subst([x=x[i-1],y=h[i-1]],rhs),
    m2:subst([x=x[i],y=h[i-1]+m1*dx],rhs),
    h[i]:h[i-1]+(m1+m2)/2*dx,

    /* Runge-Kutta */
    m1:subst([x=x[i-1],y=rk[i-1]],rhs),
    m2:subst([x=x[i-1]+dx/2,y=rk[i-1]+m1*dx/2],rhs),
    m3:subst([x=x[i-1]+dx/2,y=rk[i-1]+m2*dx/2],rhs),
    m4:subst([x=x[i],y=rk[i-1]+m3*dx],rhs),
    rk[i]:rk[i-1]+(m1+2*m2+2*m3+m4)/6*dx
)$

/* Convert exact solution and arrays to lists for graphing (starting at s) */
s:n-8;
xs:create_list(x[k],k,s,n)$
eys:create_list(subst([K=y[0],x=x[k]],es),k,s,n)$
ys:create_list(y[k],k,s,n)$
hs:create_list(h[k],k,s,n)$
rks:create_list(rk[k],k,s,n)$
plot2d([
    [discrete,xs,eys],
    [discrete,xs,ys],
    [discrete,xs,hs],
    [discrete,xs,rks]
])$