function [x,y]=rk4ode(xini,xend,yini,n,f)
%
%
%

h=(xend-xini)/n;
                                      
x = xini:h:xend;                                         
y = zeros(1,length(x)); 
y(1)=yini;
            
for i=1:(length(x)-1)                              
    k_1 = f(x(i),y(i));
    k_2 = f(x(i)+0.5*h,y(i)+0.5*h*k_1);
    k_3 = f((x(i)+0.5*h),(y(i)+0.5*h*k_2));
    k_4 = f((x(i)+h),(y(i)+k_3*h));

    y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  % main equation
end
x=x';y=y';

