function [tout, yout]= odRK4(Fcn , tspan, y0, Nstp);
t0=tspan(1);                   %tspan=[t0,tfinal] som i matlabs odelösare;
tfinal=tspan(length(tspan));
h = (tfinal - t0)/Nstp;        %steglängden beräknas
t = t0;
y = y0;                        %y0 skall vara en kolumnvektor
tout = t; 
yout = y.'                     %första raden i yout innehåller 
                               % y transponerad till radvektor.
                               % .' transponerar utan att komplex-
                               % konjugera(viktigt om y är komplex)        


for i = 1:Nstp,                %Fcn måste leverera en kolumnvektor;     
   s1=h*feval(Fcn,t,y);        %motsvarar s1=Fcn(t,y);
   s2=h*feval(Fcn,t+h/2,y+s1/2);
   s3=h*feval(Fcn,t+h/2,y+s2/2);
   s4=h*feval(Fcn,t+h,y+s3);
   y = y+(s1+2*s2+2*s3+s4)/6;  
   t = t+h;
   tout = [tout; t];           %lägg till ett värde i tout
   yout = [yout; y.'];         %lägg till en rad i yout

end;

  



