x0=0;xn=4;h=0.5; %range of solution
x=[x0:h:xn]; n=length(x);
x(1)=0; y(1)=1; %initial condition

for i=2:n
	k1=fxy(x(i-1),y(i-1));
	k2=fxy(x(i-1)+0.5*h,y(i-1)+0.5*h*k1);
	k3=fxy(x(i-1)+0.5*h,y(i-1)+0.5*h*k2);
	k4=fxy(x(i-1)+h,y(i-1)+h*k3);
	phi=(k1+2*k2+2*k3+k4)/6;
	y(i)=y(i-1)+phi*h;
end
x
y
plot(x,y)