%RK2
function [t,XM]=ex_Eulersys(fh,t0,x0,tmax,npts)
n=npts-1;
h=(tmax-t0)/n;
t(1)=t0;
XM(1,:)=x0;
for k=1:n
    t(k+1)=t(k)+h;
    x=XM(k,:);
    k1=h*fh(t(k),x);
    k2=h*fh(t(k+1),x+k1);
    XM(k+1,:)=XM(k,:)+(k1+k2)/2;
end;