%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
% Exercise 10.1.6, page 457 in 5th edition: 
%   Solve x'=(x+t)^2, with IC 
%           x(0)=-1
%   on [0,1].
% Use Taylor with derivatives up to 4th order; compare with Taylor series
% methods of orders 1,2 and3.
% Example usage: x=taylor_ex10_1_6([0 1],-1,0.01);

function [xm4,xm3,xm2,xm1,t] = taylor_ex10_1_6(int,x0,h)
a = int(1); b = int(2);
x = x0;
xm1=x; xm2=x; xm3=x; xm4 = x;
t = a;
n = round((int(2)-int(1))/h);
format long
for k = 1 : n
    xprime = (t + x).^2;
    xsecond = 2*(xprime+1).*(x+t);
    xtertius = 2*( xsecond.*(x+t) + (xprime + 1).^2);
    xquatro = 2*( xtertius.*(x+t) + 3*(xsecond.*(xprime+1))); 
    
    xm1 = xm1 + h *xprime;
    xm2 = xm2 + h *(xprime + .5*h*xsecond);
    xm3 = xm3 + h *(xprime + .5*h*(xsecond + (1/3)*h*xtertius));
    xm4 = xm4 + h *(xprime + .5*h*(xsecond + (1/3)*h*(xtertius + .25*h*xquatro)));
    t = a + k*h;
[s,y]=ode15s('exode10_1_6',[a b],x0) ;  
[ss,yy]=ode23s('exode10_1_6',[a b],x0) ;  
[sss,yyy]=ode45('exode10_1_6',[a b],x0) ;  
hold all
plot(t,xm1,'y-o');%,'LineWidth',1,'MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',1);
plot(t,xm2,'r-o')%,'LineWidth',1,'MarkerEdgeColor','k','MarkerFaceColor','b','MarkerSize',1);
plot(t,xm3,'b-*');%,'-*','LineWidth',1,'MarkerEdgeColor','k','MarkerFaceColor','b','MarkerSize',1);   
plot(t,xm4,'r-*')%,'-*','LineWidth',1,'MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize',1); 
plot(s,y)
plot(ss,yy)
plot(sss,yyy)
legend('Taylor 1','Taylor 2','Taylor 3','Taylor 4','ODE15s','ODE23s','ODE45')
end
% text(b,xm1,'\leftarrow x: 1st order Taylor',...
%      'HorizontalAlignment','left') 
% text(b,xm2,'\leftarrow x: 2nd order Taylor',...
%      'HorizontalAlignment','left') 
% text(b,xm3,'\leftarrow x: 3rd order Taylor',...
%      'HorizontalAlignment','left') 
% text(b,xm4,'\leftarrow x: 4th order Taylor',...
%      'HorizontalAlignment','left')
 %text(b,y,'\leftarrow y: <ode23> solution',...
 %   'HorizontalAlignment','left') 
fprintf('At time %1.2f, the solution of 4th order Taylor series is %1.15f\n', t,xm4)
fprintf('At time %1.2f, the solution of 3rd order Taylor series is %1.15f\n', t,xm3)
fprintf('At time %1.2f, the solution of 2nd order Taylor series is %1.15f\n', t,xm2)
fprintf('At time %1.2f, the solution of 1st order Taylor series is %1.15f\n', t,xm1)

end