Plotting the Koch Snowflake The Koch snowflake is a well known fractal. Beginning with an equilateral triangle, a smaller equilateral triangle is placed halfway along each edge of the shape. This process then repeats on each edge of the new shape. Here is my attempt to plot this in MATLAB. 

clc, clearvars, close all;

length = 1; % Original side length of triangle
theta = 60*pi/180; % Angle of the equilateral triangle in radians

% Points of the original equilateral triangle
p1 = [0,0];
p2 = [1/2*length,length*sin(theta)];
p3 = [length,0];

flake_points = [p1; p2; p3]; % Vertices of the snowflake

for iteration = 1:7
    n_points = size(flake_points,1); % No. points currently in data
    new_triangle_points = nan(n_points*3,2); % New points to be added
    for i = 1:n_points
        p = flake_points(i,:);
        if i == n_points
            vector = flake_points(1,:)-flake_points(i,:);
        else
            % Vector between next point and current point
            vector = flake_points(i+1,:)-flake_points(i,:);
        end
        b1 = p + 1/3*vector; % Base point 1 of the new triangle
        b2 = p + 2/3*vector; % Base point 2 of the new triangle
        top = p + 1/2*vector + perp(vector)*length/3*sin(theta);
        new_triangle_points(3*i-2:3*i,:) = [b1; top; b2];
    end
    flake_points(end+1:end+9,:) = nan(9,2); % Add rows to data
    new_flake_points = nan(n_points*3+n_points,2);
    n = 1;
    for i = 1:n_points
        new_flake_points(n,:) = flake_points(i,:);
        index = 3*i-2;
        new_flake_points(n+1:n+3,:) = new_triangle_points(index:index+2,:);
        n = n + 4; % Three new points have been inserted,
                    % so we now skip to the next empty spaces
    end
    flake_points = new_flake_points;
    length = 1/3*length; % Revise length
end

plot(flake_points(:,1),flake_points(:,2),'w')
hold on; % Plot last data point with 1st (ie, close the loop)
plot([flake_points(1,1),flake_points(end,1)],...
    [flake_points(1,2),flake_points(end,2)],'w')
set(gca,'color',[0.1 0.1 0.1])
axis equal;
str = sprintf('Iteration: %d',iteration);
title(str, 'FontSize', 20);


Output:



Edit - This function is also required:

function p_vec = perp(vec)
% Returns the unit vector p_vec, which is perpendicular to the input
% 2D vector vec
%
% Inputs
%
%   vec : double
%       n_dimensions x 1
%       Input vector
%
% Outputs
%   
%   p_vec : double
%       n_dimensions x 1
%       Perpendicular unit vector

p_vec(1) = -vec(2);
p_vec(2) = vec(1);

% Convert to unit vector
p_vec = p_vec / norm(p_vec);

end
Very nice script! 

I suggest you put this into a function, that takes n as an input argument instead of hardcoding it. It requires very little alterations, and will make the function a lot easier to use. You can skip all the clearing in the beginning and substitute it with:

function [] = koch_snowflake(iterations)  % or 


and add an end at the bottom of your script. 

length = 1; % Original side length of triangle


This is a no no. I strongly recommend you not to use length as a variable name as it's a very commonly used builtin function. This can cause a lot of errors that are hard to find. Note that you should always use the function numel instead of length, but it's still not a good idea to overload it.

theta = 60*pi/180; % Angle of the equilateral triangle in radians


I guess you think 60*pi/180 is easier to read than pi/3, but I disagree. If you want to use the number of degrees, and convert it to radians, you might want to use deg2rad.

I suggest substituting the iterator i to for instance ii if you ever work with complex numbers, to avoid potential errors and confusion. i is by default the imaginary unit in MATLAB (so is j, 1i, 1j). This is a matter of preference though.

Also, personal preference again, but I like it when comments are aligned (if possible):

length = 1;                   % Original side length of triangle
theta = 60*pi/180;            % Angle of the equilateral triangle in radians

% Points of the original equilateral triangle
p1 = [0,0];                   
p2 = [1/2*length,length*sin(theta)];
p3 = [length,0];              

flake_points = [p1; p2; p3];  % Vertices of the snowflake

for iteration = 1:7
    n_points = size(flake_points,1);          % No. points currently in data
    new_triangle_points = nan(n_points*3,2);  % New points to be added
...
...