function[x,t] = taylor3method(a,b,xa,f,h);

% This m-file implements a 3rd order Taylor series method for solving
%
% dx/dt = f(t,x), x(a)=xa
%
% Inputs:
%
% a  = left-hand endpoint of interval on which we approximate x
% b  = right-hand endpoint of interval on which we approximate x
% xa = initial condition
% f  = right-hand side function
% h = step size

n = floor((b-a)/h);
t = zeros(n+1,1);
x = zeros(n+1,1);
t(1) = a;
x(1) = xa;

for i = 1:n
    [xp,xpp,xppp] = feval(f,t(i),x(i));
    x(i+1) = x(i) + h*xp + 0.5*h^2*xpp + (1/6)*h^3*xppp;
    t(i+1) = t(i)+h;
end