function sol=coreuler(f,a,b,y0,N)

h=(b-a)/N;
x=a:h:b;
y=zeros(1,N+1);
y(1)=y0;
for n=1:N
    k=feval(f,x(n),y(n));
    y(n+1)=y(n)+h*feval(f,x(n)+h/2,y(n)+k*h/2);
end
sol=[x',y'];