%solve y' = f(t,y) by 2nd order Taylor method
% f(t,y) = y-t^2 +1

format long;

N = 10;
a = 0; 
b = 2; 
h= (b-a)/N;		%h is the time step.
ya = 0.5;       % initial value 
f = inline('y-t^2 + 1')      % the rate funtion  
df = inline('y-t^2+1-2*t')   % first derivative of f.  

% allocate array to save time points and solutions computated at these time
% points
T = zeros(1,N+1);
W = zeros(1,N+1);              % for saving approximate solutions
Yexact = zeros(1,N+1);         % for saving exact solutions

T(1) = a;
W(1) = ya; 
Yexact(1) = ya;

for j = 1:N, 
  tj = T(j);     
  wj = W(j);
  W(j+1) = W(j) + h*(f(tj,wj) + h/2*df(tj,wj));  %Difference Eqn for 
                                                 %2nd order Taylor method
                                                 
  T(j+1) = a + h*j;                              %advance to new time step
  Yexact(j+1) = (T(j+1)+1)^2-0.5*exp(T(j+1));    % compute exact solution
end

plot(T,W,'b--',T,Yexact,'r-',T,W,'o');
legend('Approximate','Exact');
title('Euler Approximation to dy/dt=y-t^2+1, h=0.5');
xlabel('Time');
ylabel('w(t), y(t)');

% display solutions
for j = 1:N+1, 
    disp(sprintf('t=%f,  w=%f,  y(t)=%f',T(j),W(j),Yexact(j)));
end