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

    //plot((t0:h:t0+T),z)
    
endfunction