function y = e_adams(p, h, y0, tfinal)
   N = round(tfinal/h);
   f = y = repmat(y0, 1, N+1);
   if (p == 1)
      for n = 1:N
         y(:, n+1) = y(:, n) + h*rhs(y(:, n));
      end
   end
   if (p == 2)
      f(:, 1) = rhs(y(:, 1));
      y(:, 2) = y(:, 1) + h*f(:, 1);
      for n = 2:N
         f(:, n) = rhs(y(:, n));
         y(:, n+1) = y(:, n) + h*(1.5*f(:, n) - 0.5*f(:, n-1));
      end
   end
   if (p == 3)
      f(:, 1) = rhs(y(:, 1));
      y(:, 2) = y(:, 1) + h*f(:, 1);
      f(:, 2) = rhs(y(:, 2));
      y(:, 3) = y(:, 2) + h*(1.5*f(:, 2) - 0.5*f(:, 1));
      b1 = 23./12.; b2 = -4./3.; b3 = 5./12.;
      for n = 3:N
         f(:, n) = rhs(y(:, n));
         y(:, n+1) = y(:, n) + h*(b1*f(:, n) + b2*f(:, n-1) + b3*f(:, n-2));
      end
   end
end
function yp = rhs(y)
   yp = [
         y(2);
         10*(1 - y(1)^2)*y(2) - y(1)
        ];
end