function sol=RK4(fcn,a,b,y0,N)
%generates euler soln to y'=f(x,y) in N steps
%with 4th oreder Runge Kutta algorithm
%initial condition is y(a)=y0
h=(b-a)/N;
x = a+(0:N)*h;
y(1)=y0;
for k=1:N
    k1 = feval(fcn,x(k),y(k));
    k2 = feval(fcn,x(k)+ h/2,y(k)+ h*k1/2);
    k3 = feval(fcn,x(k)+ h/2,y(k)+ h*k2/2);
    k4 = feval(fcn,x(k)+ h,y(k)+ h*k3);
    y(k+1)=y(k)+h*(k1+2*(k2 + k3) + k4)/6;
end
sol=[x',y'];