% This version solves problem 1.3.3 by Euler's method
% (blue line)
% and shows the exact solution (green line) for comparison
% To run: type 
%   [X, Y] = eul1_3(0,1,0.4,4)
%   [X, Y] = eul1_3(0,1,0.4,8)
%   [X, Y] = eul1_3(0,1,0.4,16)
%   these will plot the direction field (red lines), the exact solution
%   (green line) and the approximate solutions for the three different
%   step sizes (.1, .05, .025). 
%   Turn in the output for all cases; make a table of the values returned by the 
% three different runs, for the x values: 0.1, 0.2, 0.3, 0.4
%                ---------------------------------------
%                 x  |  y1    |   y2   |   y3   | yexact
%                ---------------------------------------
%                0.1 | 1.2500 | 1.2600 | 1.2655 | 1.2714
%                0.2 | 1.5400 | 1.5641 | 1.5775 | 1.5918
%                0.3 | 1.8780 | 1.9216 | 1.9459 | 1.9721
%                0.4 | 2.2736 | 2.3436 | 2.3829 | 2.4255
%                ---------------------------------------

function [X, Y] = eul1_3(x,y,x1,n)
  h = (x1-x)/n;
  X = x;
  Y = y;
  for i = 1:n
      k = f(x,y);
      x = x + h;
      y = y + h*k;
      X = [X; x];
      Y = [Y; y];
  end
  hold on
  plot(X,Y, 'b-')
  axis([0 0.5 1 3])
  X1 = linspace(0,0.5,21); Y1 = linspace(1,3, 21);
  hold on
%draw direction field
  df(@f,X1,Y1,'x','y','dirfield','r',0); 
  
  % plot exact solution
  yex = exact(X);
  plot(X,yex,'g-')
  % plot initial condition
  y0 = 1; x0 = 0;
  plot(x0,y0,'kp')

  % define differential equation
   function yp = f(x,y)
  % warning: f must be vectorized, ie use x.^2
    yp = 0.5 - x + 2*y;

   % define exact solution
      function ye = exact(x)
          ye = .5*x + exp(2*x);
    
