function [fx, fy] = interpHermite (p0, p1, m0, m1)
    t = linspace(1,5);
    u = (t' - t(1))./(t(end) - t(1));

    U(:, 1) = u.^3;
    U(:, 2) = u.^2;
    U(:, 3) = u;
    U(:, 4) = 1;
    A = [2 -2 1 1; -3 3 -2 -1; 0 0 1 0; 1 0 0 0];

    fx = U*A*[p0(1) p1(1) m0(1) m1(1)].';
    fy = U*A*[p0(2) p1(2) m0(2) m1(2)].';

    hold on
    scatter([p0(1) p1(1)], [p0(2) p1(2)], "filled");
    plot(fx, fy)
    hold off
end