function [X] = pointmilieu(a, t0, T, N, f)
    p = size (a, 1);
    X = zeros (p, N);
    
    h = T / N;
    t = t0;
    X(:, 1) = a;
    
    for j = 1:N-1
        K0 = f( t, X(:, j))
        K1 = f( t + h/2, X(:, j) + h * K0/2)
        X(:, j+1) = X(:, j) + h * K1
        t = t + h;
    end
    
endfunction