Plotting the Koch Snowflake
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
...
...