function [x,y]=myimpeuler(f,xinit,yinit,b,n)
     h=(b-xinit)/n;
     x=zeros(n+1,1);
     y=zeros(n+1,1);
     x(1)=xinit;
     y(1)=yinit;
     for i=1:n;
     ytemp=y(i)+h*f(x(i),y(i));
     x(i+1)=x(i)+h;
     y(i+1)=y(i)+.5*h*(f(x(i),y(i))+f(x(i+1),ytemp));
end
