function [t,u]=eestab(f,t0,T,y0,df,alpha,hmax)
t(1)=t0;
u(1)=y0;
der=df(t0,y0);
if abs(der)<=1.e-3
h=hmax;
else
h=min(hmax,2*alpha/abs(df(t0,y0)));
end
t(2)=t(1)+h;
n=1;
while t(n+1)<=T
u(n+1)=u(n)+h*f(t(n),u(n));
n=n+1;
der=df(t(n),u(n));
if abs(der)<=1.e-3
h=hmax;
else
h=min(hmax,2*alpha/abs(df(t(n),u(n))));
end
t(n+1)=t(n)+h;
end
if t(n)<T
t(n+1)=T;
h=T-t(n);
u(n+1)=u(n)+h*f(t(n),u(n));
end