function main()
  % modify these values for your problem
  n = 5;  % number of steps
  x0 = 0; % initial x-position
  y0 = 1; % initial y-value
  xf = .5; % final x-position.

  [X,Y]=simple_euler(x0,y0,xf,n);
  %disp(['X values at each step:'])
  %X
  %disp(['Y values at each step:'])
  %Y
  disp(['X, Euler Y, and exact Y values at each step:']);
  [X,Y,sol(X)]
end

% modify this function for your problem.
function out = f(x,y)
  out = x - y;
end

function out = sol(x)
  out = 2*exp(-x)+x-1;
end

function [X,Y] = simple_euler(x,y,x1,n)
  h = (x1-x)/n; % step size
  X = zeros(n,1); % array for x-positions
  Y = zeros(n,1); % array for y-values
  X(1) = x;
  Y(1) = y;
  for i=2:(n+1)
    y = y + h*f(x,y);
    x = x+h;
    X(i) = x;
    Y(i) = y;
  end
end

function [X,Y] = improved_euler(x,y,x1,n)
  h = (x1-x)/n; % step size
  X = zeros(n,1); % array for x-positions
  Y = zeros(n,1); % array for y-values
  X(1) = x;
  Y(1) = y;
  for i=2:(n+1)
    fxy = f(x,y);
    u = y + h*fxy;
    x = x+h;
    y = y + h*(fxy+f(x,u))/2.;
    X(i) = x;
    Y(i) = y;
  end
end

