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