Plot Lagrange parabolas for a Simpson's composite rule

This is my code:

clearvars; close all; clc;

f =@ (x) 2 + sin((pi/20)*x) + sin((pi/3)*x);

x0 = input("Inferior limit in x to aproximate: ") ;
xm = input("Exterior limit in x to approximate: ");
n = input("Number of parabolas for the Simpson's composite rule: ");
m = 2*n;
h = (xm-x0)/(m);

X = x0:h:xm;
Y = f(X);

Zint = zeros(1,n);
Zext = zeros(1,n+1);

evenidx =@ (v) v(2:2:end);
oddidx =@ (v) v(1:2:end);

Xint = evenidx(X);
Xext = oddidx(X);

Yint = f(Xint);
Yext = f(Xext);

fplot(f, [x0,xm], "b")      % Function f(x)
hold on
plot([Xint; Xint], [Zint; Yint], "g")   % This marks the interior lines in green
plot([Xext; Xext], [Zext; Yext], "r")   % This marks the exterior lines in red (the limits of the parabolas)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:2:m
    P =@ (x) ((x-X(i+1))*(x-X(i+2)))/((X(i)-X(i+1))*(X(i)-x(i+2)))*Y(i) + ...
        ((x-X(i))*(x-X(i+2)))/(X(i+1)-X(i))*(X(i+1)-X(i+2))*Y(i+1) + ...            % Within here lies my problem
        ((x-X(i))*(x-X(i+1)))/(X(i+2)-X(i))*(X(i+2)-X(i+1))*Y(i+2);
    flot(P, [X(i), X(i+2)], "r")
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In this code I provide the function and the intervals I want in my Simpson's rule and I expect to plot the Simpson's rule graphically. My question comes from the last point. I was trying to find out an automated way to plot the parabolas between one point X(i) and the point X(i+2) while touching the point X(i+1). For that I tried to draw an fplot of P being the Lagrange polynomial formula inside a loop. As I can see, this cannot be done in Matlab, at least in this way.

There is any way to draw the parabolas betweeen the red lines? Thank you so much!

PD: The function f(x) is totally arbitrary and can be changed as one likes

Answer 1

[Torsten]

clearvars; close all; clc;

f =@ (x) 2 + sin((pi/20)*x) + sin((pi/3)*x);

x0 = 0;%input("Inferior limit in x to aproximate: ") ;
xm = 1;%input("Exterior limit in x to approximate: ");
n = 5;%input("Number of parabolas for the Simpson's composite rule: ");
m = 2*n;
h = (xm-x0)/(m);

X = x0:h:xm;
Y = f(X);

Zint = zeros(1,n);
Zext = zeros(1,n+1);

evenidx =@ (v) v(2:2:end);
oddidx =@ (v) v(1:2:end);

Xint = evenidx(X);
Xext = oddidx(X);

Yint = f(Xint);
Yext = f(Xext);

fplot(f, [x0,xm], "b")      % Function f(x)
hold on
plot([Xint; Xint], [Zint; Yint], "g")   % This marks the interior lines in green
plot([Xext; Xext], [Zext; Yext], "r")   % This marks the exterior lines in red (the limits of the parabolas)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:2:m
    P =@ (x) ((x-X(i+1)).*(x-X(i+2)))/((X(i)-X(i+1))*(X(i)-X(i+2)))*Y(i) + ...
        ((x-X(i)).*(x-X(i+2)))/((X(i+1)-X(i))*(X(i+1)-X(i+2)))*Y(i+1) + ...            % Within here lies my problem
        ((x-X(i)).*(x-X(i+1)))/((X(i+2)-X(i))*(X(i+2)-X(i+1)))*Y(i+2);
    fplot(P, [X(i)-2, X(i+2)+2], "r")
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hold off

[Joan]

It works so fine!!!

What was the error in my code? I don't seem to find it.

Thank you so much for the help! :)

[Torsten]

Compare the definitions of P - there were several errors in yours (x instead of X, missing brackets, normal multiplication instead of elementwise multiplication).

[Joan]

I see the errors now. I guess that's what happens when I don't check number per number.

Again, thank you so much for the help!