%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.5: 
%   Solve x'=t+x+x^2, with IC 
%           x(0)=1
%   on [0,1].
% Example usage: x=taylor_ex10_1_5([0 1],1,0.01);
% For a `Final Value Problem', 
% (going backward in time from a t=1 to t=0, and final value x(1)=1)
% one should use a 'negative stepsize' and `invert' the interval, 
% e.g. taylor_ex10_1_5([1 0],1,-0.01);

function [x,y,t] = taylor_ex10_1_5(int,x0,h)
a = int(1); b = int(2);
x = x0;
t = a;
[s,y]=ode45('exode10_1_5',[a b],x0) ;   
%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 = t + x + x.^2;
    xsecond = 1+ xprime + 2*x.*xprime;
    xtertius = xsecond + 2*( xprime.^2 + x*xsecond );
    xquatro = xtertius + 2*( 3*xprime.*xsecond + x.*xtertius ); 
    xpentium = xquatro + 2*( 3*xsecond.^2 + 4*xprime*xtertius + x.*xquatro); 
    x = x + h *(xprime + .5*h*(xsecond + (1/3)*h*(xtertius + .25*h*(xquatro+h/5*xpentium))));
    t = a + k*h;  
figure(1)
hold all
title('Solution using the 5th order Taylor, ODE45 and Forward Euler')
plot(t,x,'b.-');
%                 ,'LineWidth',2,...
%                 'MarkerEdgeColor','k',...
%                 'MarkerFaceColor','g',...
%                 'MarkerSize',1);
plot(s,y,'r-');
%                ,'LineWidth',1,...
%                'MarkerEdgeColor','k',...
%                'MarkerFaceColor','r',...
%                'MarkerSize',2);
%legend('Taylor order 5','ODE45')
%end
% text(b,x,'\leftarrow x: 5th order Taylor',...
%      'HorizontalAlignment','left') 
% text(s(end),y(end),'\leftarrow y: <ode45> solution',...
%      'HorizontalAlignment','left') 

tE(1)=int(1); yE(1)=x0;

n=round((int(2)-int(1))/h);
for i=1:n
  tE(i+1)= tE(i)+h;
  yE(i+1) = eulerstep(tE(i),yE(i),h);
end
%figure(2)
hold all
plot(tE,yE,'g.-')
%title('Solution using 1st order forward Euler')
legend('Taylor order 5','ODE45','Forward Euler')
end
fprintf('At time %1.2f, the solution of 4th order Taylor series is %1.15f\n', t(1),x(end))

function y=eulerstep(t,y,h)
%one step of Euler's Method
%Input: current time t, current value y,  stepsize h
%Output: approximate solution value at time t+h
y=y+h*ydot(t,y);
end
function z=ydot(t,y)
z = t+y+y.^2;
end
end