function [t,a] = heun(h, npts, t0, y0)
    t(1) = t0;
    a(1) = y0;
    for i = 2: npts+1;
        ft   = f(t(i-1),a(i-1));
        t(i) = t(i-1) + h;
        a(i) = a(i-1) + (h/2)*(ft + f(t(i), a(i-1) + h*ft));
    end;