It looks like you are plotting your data one point at a time. MATLAB does not automatically include a marker in the line format, so when you plot a single point, you cannot see the line.
The fix here is to add a markerstyle to your plot command.
clear, clear all, clc
vacuum_permittivity = 1; %Defines epsilon 0 value = 1
for num = 1:11 %Initialises for loop to the length of number_of_point
number_of_point=11:4:51; %Initialises array of the number of points
number_of_points=number_of_point(num); %Selects the part of the number of points array to use for the conditions below
x_1 = -5; %Defines the minimum x axis value
x_N = 5; %Defines the maximum x axis value
x = transpose(linspace(x_1, x_N, number_of_points)); %Transposes x as a matrix with N number of points
y = x; %Defines the x dimensions to be equal to the y dimensions
h = x(2) - x(1); %Finds the 'h' value as the spacings between the points
[X, Y] = meshgrid(x,y); %Initialises the meshgrid for the
electricpotential_analytic=(X.^2+Y.^2).*exp(-X.^2-Y.^2); %Defines the analytical electric potential
sigma=exp(-X.^2-Y.^2).*(12*X.^2-4*X.^2.*(X.^2+Y.^2)+12*Y.^2-4*Y.^2.*(X.^2+Y.^2)-4); %Defines the analytical charge density
% Define our differential operator for two dimensions as a matrix.
A = (1 / h^2) * (diag(ones(1, number_of_points - 1), 1) + ...
diag(-4 * ones(1, number_of_points), 0) + ...
diag(ones(1, number_of_points - 1), -1));
B = (1 / h^2) * diag(ones(1, number_of_points), 0);
L = kron(diag(ones(1, number_of_points), 0), A) + ...
kron(diag(ones(1, number_of_points - 1), 1), B) + ...
kron(diag(ones(1, number_of_points - 1), -1), B);
sigma = reshape(sigma, number_of_points^2, 1); %Reshapes sigma to have the same dimensions as the electric potential
electricpotential_analytic=reshape(electricpotential_analytic,number_of_points^2,1); %Reshapes the electric potential to have the same dimensions as sigma for a matrix
electric_potential=-inv(L)*sigma; %Uses equation 39 to define Poisson's equation in matrix format. No element wise multiplication is done since are interested in the error in scalar format and not as a matrix
error=max(abs(electricpotential_analytic-electric_potential))/(max(abs(electricpotential_analytic))) %Determines the error from equation 50 as a scalar value
% Add markerstyle -------> vvv
plot(number_of_points,error,'o') %Plots the error function for every number of points
hold on %Holds the value so that multiple values can be plotted on the same figure
end %Terminates the for loop
xlabel('N','fontsize', 16) %Labels the x axis
ylabel('Error','fontsize', 16) %Labels the y axis
