x0=0;xn=4;h=0.5; %range of solution
x=[x0:h:xn]; n=length(x);
x(1)=0; y(1)=2; %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+1/sqrt(2))*h*k1+(1-1/sqrt(2))*h*k2);
	k4=fxy(x(i-1)+h,y(i-1)-1/sqrt(2)*h*k2+(1+1/sqrt(2))*h*k3);
	phi=(k1+2*(1-1/sqrt(2))*k2+2*(1+1/sqrt(2))*k3+k4)/6;
	y(i)=y(i-1)+phi*h;
end
x;
y;
plot(x,y)