function [t,w]=rk4(y0,a,b,h)
% rk4(1,0,1,0.1);

N=(b-a)/h;
t=linspace(a,b,N+1);
w(1)=y0;
for i=1:N
    k1=h*defunc(t(i),w(i));
    k2=h*defunc(t(i)+h/2,w(i)+k1/2);
    k3=h*defunc(t(i)+h/2,w(i)+k2/2);
    k4=h*defunc(t(i)+h,w(i)+k3);
    w(i+1)=w(i)+(k1+2*k2+2*k3+k4)/6;
end
[T,Y] = ODE45('defunc',t,y0);
plot(T,Y)
hold on
plot(t,w,'o')
hold off
error=max(abs(w-Y'))
