function [t,y] = rk(f,tint,y0,tol,h0)
% 
%  [t,y] = rk(f,tint,y0,tol,h0)
%  
%  Loeser y'=f(t,y) med en et 1,2 Runge-Kutta par. 
%  Input: f   : Haandtak til funksjonen f
%         tint: Integrasjonsintervallet
%         y0  : Startverdier
%         tol : Toleranse
%         h0  : Skrittlengde for foerste skritt. 

t0 = tint(1);
tend = tint(2);

% Arbeidsvariable
tn = t0;
yn = y0; 
h = h0;

% Lagrer loesningne for hvert skritt i t og y. 
t = t0;
y = y0;

while (tn<tend)
    if (tend-tn < h)
        h = tend-tn;
    end
    k1 = feval(f,tn,yn);
    k2 = feval(f,tn+h,yn+h*k1);
    feil = norm(0.5*h*(-k1+k2));
    if feil<tol
        % Skrittet aksepteres. 
        tn = tn+h;
        yn = yn+0.5*h*(k1+k2); % Lokal ekstrapolasjon
        
        y = [y,yn];            % Lagrer loesningen
        t = [t,tn];
    end
    % Skrittlengden forandres. 
    fac = 0.75*sqrt(tol/feil);
    fac = min(5,max(0.1,fac)); 
    h = fac*h;
end

    
    
