function [X,Y]=impeuler(x,y,x1,n)

h=(x1-x)/n;         % step size
X=x;                % initial x
Y=y;                % initial y

for i=1:n           % begin loop
   k1=f(x,y);       % first slope
   k2=f(x+h,y+h*k1);% second slope
   k=(k1+k2)/2;     % average slope
   x=x+h;           % new x
   y=y+h*k;         % new y
   X=[X,x];         % update x-column
   Y=[Y,y];         % update y-column
end                 % end loop