What is the best way to calculate and store coordinate data for a multi-link mechanism?

I want to calculate the coordinates (x,y) of a number of points inside the working area of a mechanism in order to later determine the maximum load the mechanism can lift with its endpoint.

I have the basic equations figured out for a single point using position vectors and translation matrices but am having problem find a way to conveniently store and then plot the data to visually verify my results.

For my basic 2-link test I used a double for loop and stored the x and y coordinates in separate mxn arrays (where m is the range of motion for the first link and n the range of motion for the second link, or the other way around... irrelevant). Which works but isn't a very elegant solution.

The problem comes when adding a third link. lugging around with a bunch of 3D arrays quickly becomes cumbersome. (for the final program I also need to know the positions of a few CG's and some other interesting points) Anyone have any suggestions on how to proceed?

The solution is probably hilariously simple it's just that I'm too blind, green and tired to see it.

Above: working area for 2-link test.

Below: Relevant piece of code.

%Geometry model
i=1;
k=1;
  %Pre-allocating matrices
  xTIP = zeros(Steps2+1,Steps1+1);
  yTIP = zeros(Steps2+1,Steps1+1);

 for CL1 = CL1_min:Step_size1:CL1_max;
    phi2 = acos((norm(rC1PJ1)^2+norm(rLoC1L)^2-CL1^2)/(2*norm(rC1PJ1)*norm(rLoC1L)));
    thetaL = phi1 + phi2 + phi3; %Angle between reference systems of link 1 (P) and 2 (L)

    APL = [cos(pi/2-thetaL) sin(pi/2-thetaL); -sin(pi/2-thetaL) cos(pi/2-thetaL)]; %Transformation matrix, rotation around J1/Lo

    for CL2 = CL2_min:Step_size2:CL2_max;
        phi4 = acos(((norm(rC2LJ2)^2)+(norm(rVoC2V)^2)-(CL2^2))/(2*norm(rC2LJ2)*norm(rVoC2V)));
        thetaV = phi4 + phi5 + phi6; %Angle between reference systems of link 2 (L) and link 3 (V)

        ALV = [cos(pi-thetaV) sin(pi-thetaV);-sin(pi-thetaV) cos(pi-thetaV)]; %Transformation matrix, rotation around J2/Vo

        rTIP = rPoJ1+APL*rLoJ2+APL*ALV*rVoJ3; %Position vector for tip with respect to the mechanism origin.
        xTIP(k,i) = rTIP(1);
        yTIP(k,i) = rTIP(2);

        k=k+1;
    end
    k=1;
    i=i+1;
end

 %Plot working area
figure(1)

    %Plot positions
    plot(xTIP(2:end-1,2:end-1),yTIP(2:end-1,2:end-1),'o');
    hold on
    %Plot working area boundaries
    plot(xTIP(1,:),yTIP(1,:),'black',xTIP(end,:),yTIP(end,:),'black',xTIP(:,1),yTIP(:,1),'black',xTIP(:,end),yTIP(:,end),'black')