% To run: type [X, Y] = impeuler(x0,y0,xfinal,steps)
% enter this command repeatedly for different
% x0, y0 to get plots of multiple solutions on the 
% same graph. Or, you may keep x0, y0 fixed and vary
% the number of steps to get a feeling of what 
% is the error involed in the computations.
% To change the problem solved, simply edit the
% function f(x,y)

function [X, Y] = impeuler(x,y,x1,n)
  h = (x1-x)/n;
  X = x;
  Y = y;
  for i = 1:n
      k1 = f(x,y);
      k2 = f(x+h, y+h*k1);
      k = (k1 + k2)/2;
      x = x + h;
      y = y + h*k;
      X = [X; x];
      Y = [Y; y];
  end
  hold on
  plot(X,Y)
  axis([-2 2 -2 2])
  X1 = linspace(-2,2,21); Y1 = linspace(-2,2,21);
  hold on
  df(@f,X1,Y1,'x','y','dirfield','r',0); 
  
  function yp = f(x,y)
    yp = 3-y;
    
    
