n=(tf-t0)/h;
t=t0;
x=x0;
y=y0;
tvect=[t0];
xvect=[x0];
yvect=[y0];
for i=1:n
   xnew=eval(eulerfcn1);
   ynew=eval(eulerfcn2);
   t=t+h;
   x=xnew;
   y=ynew;
   tvect=[tvect,t];
   xvect=[xvect,x];
   yvect=[yvect,y];
end
tvect
xvect
yvect
plot(tvect,xvect,tvect,yvect)
xlabel('t')
ylabel('x,y')