function [t,y]=fe(f,tspan,y0,k)
%Example:
%f=@(t,y) -5*t.*y.^2+5./t-1./(t.^2);
%[t,y]=fe(f,[1 25],1,.21);
%plot(t,y)
%[t,y]=fe(f,[1 25],1,.19);
%hold on
%plot(t,y)
  
  
  T=tspan(2)-tspan(1);
  t=[tspan(1):k:tspan(2)];
  N=length(t);
  y=zeros(N,1);
  y(1)=y0;
  for i=2:N,
    y(i)=y(i-1)+k*feval(f,t(i-1),y(i-1));
  end
  
