function [t, u] = rk2(f, tspan, u0, n)
% RK2    Second-order Runge-Kutta for an IVP.
% Input:
%   f       defines f in u'(t)=f(t,u) (function)
%   tspan   endpoints of time interval (2-vector)
%   u0      initial value (vector, length m)
%   n       number of time steps (integer)
% Output:
%   t       selected nodes (vector, length n+1)
%   u       solution values (array, (n+1)-by-m)

    % Time discretization.
    a = tspan(1);  
    b = tspan(2);
    h = (b-a)/n;   
    t = a + h*(0:n)';
    
    % Initialize solution array. 
    u = zeros(length(u0), n+1);  
    u(:,1) = u0;
    
    % Time stepping. 
    for i = 1:n 
        k1 = h*f( t(i),     u(:,i) );
        k2 = h*f( t(i)+h,   u(:,i)+k1 );
        u(:,i+1) = u(:,i) + (k1 + k2)/2;
    end

    u = u.';   % conform with MATLAB output convention
end