%Taylor series Method of higher order for Solving Initial Value Problems
% Input: interval [a,b], initial value x0, step size h
% Output: time steps t, solution y
% Example usage: x=taylor([1 2],2,0.01);

function [x,y,t] = taylor(int,x0,h)
a = int(1); b = int(2);
x = x0;
y = x0;
t = a;
plot(t,x,'*','LineWidth',6);
plot(t,y,'o','LineWidth',3);
n = round((int(2)-int(1))/h);
format long
for k = 1 : n
    xprime = 1 + x.^2 + t^3;
    yprime = 1 + y.^2 + t^3;
    xsecond = 2*x.*xprime + 3*t^2;
    xtertius  = 2*x*xsecond + 2* xprime*xprime + 6*t;
    xquatro = 2*x*xtertius + 6* xprime*xsecond + 6;
    x = x + h *(xprime + .5*h*(xsecond + (1/3)*h*(xtertius + .25*h*xquatro)));
    y = y + h*yprime;   % Euler Method
    t = a + k*h;
hold all
plot(t,x,'*','LineWidth',3,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',10);
plot(t,y,'o','LineWidth',1,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',10);
end
text(b,x,'\leftarrow x: 4th order Taylor',...
     'HorizontalAlignment','left') 
text(b,y,'\leftarrow y: 1st order Taylor',...
     'HorizontalAlignment','left') 