function [T,Y] = euleri(f,a,b,ya,m)
% [T,Y] = euleri(f,a,b,ya,m)
% Implicit Euler's solution for y' = f(t,y) with y(a) = ya.
% f  is the function, input.
% a  is the left endpoint, input.
% b  is the right endpoint, input.
% ya is the initial condition, input.
% m  is the number of steps, input.
% T  is the vector of abscissas, output.
% Y  is the vector of ordinates, output.
h = (b - a)/m;
T = zeros(1,m+1);
Y = zeros(1,m+1);
T(1) = a;
Y(1) = ya;
for j=1:m,
  yo = Y(j);
  fo = h*feval(f,T(j),yo);
  y1 = yo + fo;
  tj1 = a + h*j;
  for k=1:5
     f1 = h * feval(f,tj1,y1);
     y1 = y1 - (y1-yo) * (y1 - Y(j) - f1) / ((y1-yo) - (f1-fo));
     if (y1 == yo), break;
     end
     yo = y1;
     fo = f1;
  end
  Y(j+1) = y1;
  T(j+1) = tj1;
end
