function [t,x] = taylor4(tRange, x0, n)
% taylor4: 4th order Taylor method for solving a single ODE
% tRange: [t0,tEnd], the starting and ending times
% x0: initial value of x
% n: number of time steps so time values are t(1) to t(n+1)
% [t,x]: an m by 2 matrix containing the (t,x) values
% t and x are each column vectors.

t0 = tRange(1);
tEnd = tRange(2);
h = (tEnd - t0) / n;
t = (t0 : h : tEnd)'; % column vector
x = zeros(length(t),1);

x(1) = x0;
for k = 1 : length(t)-1
   x1 = 1 + x(k)^2 + t(k)^3;
   x2 = 2*x(k)*x1 + 3*t(k)^2;
   x3 = 2*x(k)*x2 + 2*x1*x1 + 6*t(k);
   x4 = 2*x(k)*x3 + 6*x1*x2 + 6;
   x(k+1) = x(k) + h*(x1 + (h/2)*(x2 + (h/3)*(x3 + (h/4)*x4)));
end