function RK_test(f,ex,y0,t0,t1,N)
cla;
clf;
%f=@(t,y)(2*pi*cos(2*pi*t)); % RHS
% N=9; % number of steps
% t0=0; % start
% t1=3; % stop
h=(t1-t0)/N;
% y0=1;

t=linspace(t0,t1,N+1);
%tt=t;
tt=linspace(t0,t1,10*N+10);

%ex=@(t)(1./(t+1/y0));
%ex=@(t)(0.5*t.^2+1);
%ex=@(t,t0,y0)(1./(t-t0+1/y0));
%ex=@(t)(y0*exp(t));
%ex=@(t)(sin(t+asin(y0)));
%ex=@(t)(1+sin(2*pi*t));
%ex=@(t)(exp(sin(2*pi*t)));

% numerical solution
y=zeros(N+1,1);
y(1)=y0;

% exact solution
ey=ex(tt);

% semi-exact solution
sey=zeros(N+1,1);
sey(1)=y0;

%axis([t0 t1 0 1]);
hold on;
for l=1:N
    yl=y(l);
    tl=t0+(l-1)*h;
    
    % numerical k_i's
    k_1=f(tl,yl);
    k_2=f(tl+0.5*h,yl+0.5*h*k_1);
    k_3=f(tl+0.5*h,yl+0.5*h*k_2);
    k_4=f(tl+h,yl+h*k_3);
        
    y(l+1)=yl+h*(k_1+2*k_2+2*k_3+k_4)/6;
    
    %line([tl tl+h],[yl y(l+1)],'Color','m');
    plot(tl,yl,'*m');
    
    
    % semi-exact k_i's
    seyl=sey(l);
    sek_1=f(tl,seyl);
    sek_2=f(tl+0.5*h,ex(tl+0.5*h));
    sek_3=f(tl+0.5*h,ex(tl+0.5*h));
    sek_4=f(tl+h,ex(tl+h));
        
    sey(l+1)=seyl+h*(sek_1+2*sek_2+2*sek_3+sek_4)/6;
    
    %line([tl tl+h],[seyl sey(l+1)],'Color','k');
    plot(tl,seyl,'ok');
    
    line(tt((10*(l-1)+1):(10*l+1)),ey((10*(l-1)+1):(10*l+1)),'Color','g');
    
    %drawnow;
    %pause(0.5);
    
    
    
    %local_error=eyl1-ex(tl+h);
    %fprintf('delta k1: %0.3g\ndelta k2: %0.3g\ndelta k3: %0.3g\ndelta k4: %0.3g\n',(k_1-sek_1)/sek_1,(k_2-sek_2)/sek_2,(k_3-sek_3)/sek_3,(k_4-se_4)/se_4);
    %fprintf('Local Error: %0.3g (time %d)\n',local_error,l);
end
plot(t(end),y(end),'*m');
plot(t(end),sey(end),'ok');
line(tt((10*N+1):(10*N+10)),ey((10*N+1):(10*N+10)),'Color','g');
legend({'Numerical (Runge Kutta)','Semi-exact','Exact'});

hold off;


%plot(tt,ex(tt),t,sey,t,y);
