Continious Bezier Curves in Matlab (iterative)

%% ------------------------------------------------------------------------
% Cleanup
%%
clc;
clear variables;
figure, hold on;
%axis([0 10 -1 9]);

%% ------------------------------------------------------------------------
% Define
%%
n = 4;
pmts = zeros([n 5]);

pnts(1, 1:2) = [1 0];   % Point (x|y)
pnts(1, 3)   = 7.8*pi/8; % Angle
pnts(1, 4:5) = [0 3];   % Vector lengths

pnts(2, 1:2) = [9 1];  % Point (x|y)
pnts(2, 3)   = -7*pi/8; % Angle
pnts(2, 4:5) = [3 1.5];  % Vector lengths

pnts(3, 1:2) = [12 5];  % Point (x|y)
pnts(3, 3)   = -pi; % Angle
pnts(3, 4:5) = [3 1.5];  % Vector lengths

pnts(4, 1:2) = [17 3];  % Point (x|y)
pnts(4, 3)   = 5*pi/8; % Angle
pnts(4, 4:5) = [1.5 0];  % Vector lengths

%% ------------------------------------------------------------------------
% Draw
%%
plot(pnts(:, 1), pnts(:, 2), 's:', 'Color', [1 0 0]);
for i=1:n
    text(pnts(i, 1)+0.2, pnts(i, 2)+0.2, sprintf('P%d', i));
    for j=4:5
        if pnts(i, j) > 0
            angle = pnts(i, 3);
            if j==5
                angle = angle + pi;
            end
            length = pnts(i, j);
            ctrl_x = pnts(i, 1)+cos(angle)*length;
            ctrl_y = pnts(i, 2)+sin(angle)*length;
            plot(ctrl_x, ctrl_y, 'o', 'Color', [0.6 0.8 0.3]);
            plot([pnts(i, 1) ctrl_x], [pnts(i, 2) ctrl_y], '-', 'Color', [0.6 0.8 0.3]);
            text(ctrl_x+0.2, ctrl_y+0.2, sprintf('P%d_%d', i, j-3));
        end
    end
end
%% ------------------------------------------------------------------------
% !)(&/§%!"&/$§/&"%$
%%

% number of curves
N=n-1;
% for each ancient two points (draw a bezier curve)
for q=1:N
    %result vector
    B_x=[];
    B_y=[];
    % ttl points (pnts + ctrlpnts)
    z = 2;
    k = (z-2)*2 + 2 + z;
    % point buffer
    D = zeros(k, k, 2); % step * point * (x,y)
    % initialize from points, Pa_0, Pa, Pa_1, Pb_0, Pb, Pb_1 etc...

    % pnt q
    D(1, 1, :) = pnts(q, 1:2);
    % ctrlpng after pnt q
    angle = pnts(q, 3) + pi; % 180° ahead of ctrl pnt before pnt
    length = pnts(q, 5);
    ctrl_x = pnts(q, 1)+cos(angle)*length;
    ctrl_y = pnts(q, 2)+sin(angle)*length;
    D(1, 2, :) = [ctrl_x ctrl_y];

    % ctrlpnt before pnt q+1
    angle = pnts(q+1, 3);
    length = pnts(q+1, 4);
    ctrl_x = pnts(q+1, 1)+cos(angle)*length;
    ctrl_y = pnts(q+1, 2)+sin(angle)*length;
    D(1, 3, :) = [ctrl_x ctrl_y];

    % pnt q+1
    D(1, 4, :) = pnts(q+1, 1:2);

    % t from 0 to 1
    for t=0:0.01:1
        % One calculation step for each additional point (k-1 steps,
        % k=pnts+ctrlpnts); 1->2, 2->3, 3->4
        for j=2:k
            % calculate each line in that iteration (as many lines to calculate as iterations left)
            for i=1:k-j+1
                a = D(j-1,i,:);   % from pnt (previous iteration)
                b = D(j-1,i+1,:); % to pnt (previous iteration)
                % calc point on line
                v = a * (1-t) + b * t;
                D(j,i,:)=v;       % save for next iteration
            end
        end
        % D(k-1,1) is point on curve
        B_x=[B_x D(k,1,1)]; % X
        B_y=[B_y D(k,1,2)]; % Y
    end

    plot(B_x, B_y, 'b');
end