function [Yhat,Yn]=rk45(x,Y,h)
k1=h*f(x,Y);
k2=h*f(x+.25*h,Y+.25*k1);
k3=h*f(x+3*h/8.,Y+3*k1/32.+9*k2/32);
k4=h*f(x+12*h/13.,Y+1932*k1/2197.-7200*k2/2197.+7296*k3/2197);
k5=h*f(x+h,Y+439*k1/216.-8*k2+3680*k3/513.-845*k4/4104);
k6=h*f(x+.5*h,Y-8*k1/27.+2*k2-3544*k3/2565+1859*k4/2565-11*k5/40.);	
Yhat=Y+25*k1/216.+1408*k3/2565.+2197*k4/4104.-k5/5;
Yn=Y+16*k1/135.+6656*k3/12825+28561*k4/56430-9*k5/50+2*k6/55;