How to define NACA 5 digit series equation?

I've been trying to analytically define the shape of the NACA 5 digit series (both reflex and non-reflex). However, when trying to produce the plots, the shape that MATLAB produces is very strange. For some airfoils, the shape it produces is perfect, but for others it is not. I have included images of both. I have yet to find a pattern to see which ones fail and which don't, as well as a reason why.

Geometric definition function code:

function [x_upper, y_upper, x_lower, y_lower] = NACAFoil5(NACA5, Chord) %definition of function
x = 0:0.00001:Chord; % definition of x coordinates
t = str2num([num2str(NACA5(4)),num2str(NACA5(5))])/100; %#ok<ST2NM>
y_t = 5*t*Chord*(0.2969*sqrt(x/Chord)-0.126*(x/Chord)-0.3516*(x/Chord).^2+0.2843*(x/Chord).^3-0.1015*(x/Chord).^4); 
p = NACA5(2)/20;

if length (NACA5) == 5
    disp(['NACA 5 Digit Series: NACA ', num2str(NACA5(1)) num2str(NACA5(2)) num2str(NACA5(3)) num2str(NACA5(4)) num2str(NACA5(5))])
    
    if NACA5(3) == 0
        reflex = 0;
    elseif NACA5(3) == 1
        reflex = 1;
    end
end

if reflex == 0
    if NACA5(2) == 1
            disp('5% non reflexed camber airfoil')
            r1 = 0.0580;
            k1 = 361.4;
        
        elseif NACA5(2) == 2
            disp('10% non reflexed camber airfoil')
            r1 = 0.126;
            k1 = 51.64;
        
        elseif NACA5(2) == 3
            disp('15% non reflexed camber airfoil')
            r1 = 0.2025;
            k1 =  15.957;
        
        elseif NACA5(2) == 4
            disp('20% non reflexed camber airfoil')
            r1 = 0.29;
            k1 = 6.643;
        
        elseif NACA5(2) == 5
            disp('25% non reflexed camber airfoil')
            r1 = 0.391;
            k1 = 3.23;
     end
         for i = 1:length(x)
             if (x(i)/Chord) <= p
                 y_camber(i) = (k1/6)*(((x(i)/Chord)^3)-(3*r1*(x(i)/Chord)^2)+(r1^2)*(3-r1)*(x(i)/Chord));
                 dy_camber(i) = (k1/6)*((3*(x(i)/Chord)^2)-(6*r1*(x(i)/Chord))+(r1^2)*(3-r1));
             elseif (x(i)/Chord) > p
                 y_camber(i) = ((k1*r1^3)/6)*(1-(x(i)/Chord));
                 dy_camber(i) = (-1*k1*r1^3)/6;
             end
         end


elseif reflex == 1
    if NACA5(2) == 2
        disp('10% reflexed camber airfoil')
        r1 = 0.13;
        k1 = 51.99;
        k2_k1 = 0.000764;
        
    elseif NACA5(2) == 3
        disp('15% reflexed camber airfoil')
        r1 = 0.217;
        k1 = 15.793;
        k2_k1 = 0.00677;
        
    elseif NACA5(2) == 4
        disp('20% reflexed camber airfoil')
        r1 = 0.318;
        k1 = 6.52;
        k2_k1 = 0.303;
        
    elseif NACA5(2) == 5
        disp('25% reflexed camber airfoil')
        r1 = 0.441;
        k1 = 3.191;
        k2_k1 = 0.1355;
    end
    
    for i = 1:length(x)
        if (x(i)/Chord) <= p
            y_camber(i) = (k1/6)*((((x(i)/Chord)-r1)^3)-(k2_k1)*((1-r1)^3)*x(i)-(r1^3)*((x(i)/Chord))+r1^3);
            dy_camber(i) = (k1/6)*(3*(((x(i)/Chord)-r1)^2)-(k2_k1)*((1-r1)^3)-r1^3);
        elseif (x(i)/Chord) > p
            y_camber(i) = (k1/6)*((k2_k1)*(((x(i)/Chord)-r1)^3)-(k2_k1)*((1-r1)^3)*(x(i)/Chord)-(r1^3)*((x(i)/Chord)+r1^3));
            dy_camber(i) = (k1/6)*((3*k2_k1)*(((x(i)/Chord)-r1)^2)-(k2_k1)*((1-r1)^3)-r1^3);
        end
    end     
end
       
theta = atan(dy_camber);
         x_upper = x - (y_t.*sin(theta));
         y_upper = y_camber + (y_t.*cos(theta));
         x_lower = x + (y_t.*sin(theta));
         y_lower = y_camber - (y_t.*cos(theta));  
         
         y_upper(end) = 0;           
         y_lower(end) = 0; 
end

Solver:

clc
clear all
close all


%% Body Equation
AirfoilTypeAsk = 'Enter "J" for Jokowski Airfoil, "N5" for NACA 5 series airfoils, "N" for NACA Airfoils or "C" for Cylinder:  ';      % Body Type
Type = input(AirfoilTypeAsk);
if Type == 'N5'
AirfoilAsk = 'Enter The Airfoil Number in Raw Vector as [x x x x x]:  ';                                 % Airfoil number
Airfoil = input(AirfoilAsk);
NACA5 = Airfoil;
ChordAsk = 'Enter The Airfoil Chord Length:  ';                                                        % Airfoil Chord
Chord = input(ChordAsk);
[x_upper,y_upper,x_lower,y_lower] = NACAFoil5(NACA5,Chord);
elseif Type == 'N'
AirfoilAsk = 'Enter The Airfoil Number in Raw Vector as [x x x x]:  ';                                 % Airfoil number
Airfoil = input(AirfoilAsk);
NACA = Airfoil;
ChordAsk = 'Enter The Airfoil Chord Length:  ';                                                        % Airfoil Chord
Chord = input(ChordAsk);
[x_upper,y_upper,x_lower,y_lower] = NACAFoil(NACA,Chord);
elseif Type == 'J'
AirfoilAsk = 'Enter The Airfoil max. Thickness/Chord and max. Camber/Chord as [t/c C/c]:  ';           % Airfoil Thickness and Camber
Airfoil = input(AirfoilAsk);
ChordAsk = 'Enter The Airfoil Chord Length:  ';                                                        % Airfoil Chord
Chord = input(ChordAsk);
t_c = Airfoil(1);
C_c = Airfoil(2);
[x_upper,y_upper,x_lower,y_lower] = JFoil(t_c,C_c,Chord);    
elseif Type == 'C'
CylinderAsk = 'Enter The Radius:  '; 
Raduis = input(CylinderAsk);
Chord = 2*Raduis;
[x_upper,y_upper,x_lower,y_lower] = cylinder(Chord);    
end


figure      % Airfoil
hold on
grid on
axis equal
plot(0:0.0001:Chord,polyval(polyfit(x_upper,y_upper,10),0:0.0001:Chord),0:0.0001:Chord,polyval(polyfit(x_lower,y_lower,10),0:0.0001:Chord),'LineWidth',1.5,'color','b')
plot(0:0.0001:Chord,(polyval(polyfit(x_upper,y_upper,10),0:0.0001:Chord)+polyval(polyfit(x_lower,y_lower,10),0:0.0001:Chord))/2,'--','LineWidth',1.5,'color','r')
if Type == 'N'
title(['NACA ',  num2str(NACA(1)) num2str(NACA(2)) num2str(NACA(3)) num2str(NACA(4)),'  Airfoil'])
elseif Type == 'J'
title(['Flow Stream Lines around  Joukowski Airfoil:    t_{max}/c = ',  num2str(t_c*100),'%    C_{max}/c = ' ,num2str(C_c*100),'%']) 
elseif Type == 'N5'
title(['NACA5 ', num2str(NACA5(1)) num2str(NACA5(2)) num2str(NACA5(3)) num2str(NACA5(4)) num2str(NACA5(5)), ' Airfoil'])
end
xlabel('x')
ylabel('y')
legend('Upper Surface','Lower Surface','Camber Line')

The data, including the equations used are from the following references:

http://airfoiltools.com/airfoil/naca5digit

https://people.clarkson.edu/~pmarzocc/AE429/The%20NACA%20airfoil%20series.pdf

Also, the main code is designed to solve the flow over an arbitrary shape, I'm just trying to add the shape of NACA 5 into the solver. So kindly ignore the other shape types and their respective commands.