function [t,y] = EulersysV2(xpsys,t0,y0,tmax,npts)
n=npts-1;
h=(tmax-t0)/n;
t(1)=t0;
y(1,:)=y0;
for k=1:n
    x=y(k,:);
    f1=h*xpsys(t(k),x);
    y(k+1,:)=y(k,:)+f1;
    t(k+1)=t(k)+h;
end