The issue is x, f(x), and y are all discrete points every loop. They have no vectorized connection to their previous point. One solution I can see is only plot the result afterwards and have the x, f(x), and y be vectorized.
dy = @(x,y)y-x^2+1; % differential equation
f = @(x)((x+1)^2)-0.5*exp(x); % analytical solution
x0 = 0; % initial point of interval
xf = 2; % final point of interval
y = 0.5; % initial condition of value of y at x0
h = 0.1; % step size for side edge
fprintf('x(i)\t\t y_Euler(i)\t\t y_analitical(i)\n') % data table header
fprintf('%f\t %f\t\t %f\n',x0,y,f(x0));
count=1;
for x = x0 : h : xf-h
y = y + dy(x,y)*h; % Euler formula
x = x + h; % here we compute x as the next step
fprintf('%f\t %f\t\t %f\n',x,y,f(x)); % print results
xplot(count)=x;
yplot(count)=y;
funplot(count)=f(x);
count=count+1;
end
figure(1)
plot(xplot,funplot,'rx-','LineWidth',2)
hold on
plot(xplot,yplot,'bo')
title('Euler`s method and exact solution')
xlabel('x','FontSize',12,'FontWeight','bold','Color','black')
ylabel('y','FontSize',12,'FontWeight','bold','Color','black')
legend('Analytical solution','Euler method')
grid on; grid minor
