function z = eulertaylor(N);
% EULERTAYLOR   Plot and compare Euler's method and Taylor
%   order 2 and improved Euler on the problem
%       y' = x - e^y,   y(0) = 1.
%   Approximates y(2).  Uses N steps.
%   Example:  >> eulertaylor(10)

h = 2 / N;
 
x = 0:h:2; % equally-space x values on interval [0,2]

% space for the approximate solutions
yE = zeros(size(x));  yT = yE;  yIE = yE;

yE(1) = 1;  yT(1) = 1; yIE(1) = 1; % initial condition

for j=1:N
  % Euler:
  fE   = x(j) - exp(yE(j));    % = y' at (x_n,y_n)
  yE(j+1) = yE(j) + h * fE;

  % Taylor order 2:
  fT   = x(j) - exp(yT(j));    % = y' at (x_n,y_n)
  yppT = 1 - exp(yT(j)) * fT;  % = y'' at (x_n,y_n)
  yT(j+1) = yT(j) + h * fT + (h*h/2) * yppT;

  % improved Euler:
  fIE     = x(j) - exp(yIE(j));  % = y' at (x_n,y_n)
  ynext   = yIE(j) + h * fIE;
  fnextIE = x(j+1) - exp(ynext); % = y' at Euler step location
  yIE(j+1) = yIE(j) + h * 0.5 * (fIE + fnextIE);
end

plot(x,yE,'-*',x,yT,'-o',x,yIE,'-s')
legend('Euler','Taylor order 2','improved Euler')

ygood = 0.371389615;   % estimated by "eulertaylor(40000)"
hold on, plot(2,ygood,'ko',"markersize", 4)
text(1.5,ygood,"exact y(2) in circle"), hold off

z = [yE(N+1); yT(N+1); yIE(N+1)];