function y0=Runge_kutta(f,h,y0,x0)

syms x y

k1 = h * subs(f,{x,y},{x0,y0});
k2 = h * subs(f,{x,y},{x0+h/2,y0+k1/2});
k3 = h * subs(f,{x,y},{x0+h/2,y0+k2/2});
k4 = h * subs(f,{x,y},{x0+h,y0+k3});
y0 = y0 + 1/6*(k1+2*k2+2*k3+k4);