Maximum modulus of deviation of the value of the interpolation polynomial from the value of the original function

I am learning computations on GNU octave and am struggling with the following task.

Let the interpolation nodes X be given and the functions lagrange(X,Y,x0) and original_function(x0) defined, which calculate the values of the interpolation polynomial and the original interpolated function at the point x0. And also the delta_max(Y1,Y2) function is already defined.

I need to write code that looks at each point in the interval [xmin,xmax] with step h=xmax−xmin / 100N (N is the number of interpolation intervals) and finds the one where the maximum modulus of the deviation of the value of the interpolation polynomial from the value of the original function is reached, as well as the value of this maximum itself.

The input/output format for the task:

Sample Input:

1 2 3

Sample Output:

h = 0.010000
x_delta_max = 2.6100
delta = 0.58020

I am struggling to provide the answer in the given time (13 seconds). Is there something obvious I am missing?

I need to do this task using GNU Octave, but if there are any suggestions regarding overall structure of the code (e. g. the mistake I make is not octave-specific), MatLab examples will be useful.

First, I calculated the step and an array of all the points from this interval. x_h is for X coords, y_h is for Y. In the next step I tried to calculate the value of interpolation polynomial in every given point (as stated in the task). I think i am doing something wrong in this very moment, as the last operation (delta_max) is relatively simple.

h = (X(end) - X(1)) / (100 * (size(X, 2) - 1));
x_h = X(1):h:X(end);
y_h = arrayfun(@original_function, x_h);
lagr = arrayfun(@(x) lagrange(x_h, y_h, x), x_h);
[delta, x_delta_max] = delta_max(y_h, lagr);

But this code takes too long to execute and I am always out of any time limitations.

Below is the code for lagrange(X, Y, x0), original_function(x0) and delta_max(Y1, Y2).

lagrange() function computes the value of an interpolation polynomial at a given point. X - x coordinates of interpolation nodes, Y - y coordinates of interpolation nodes, x0 - point where the value is calculated

function y0 = lagrange(X,Y,x0)
    y0 = 0;
    n = size(X, 2);
    for i = 1:n
        p = 1;
        for j = 1:n
            if (j != i)
                p *= (x0 - X(j)) / (X(i) - X(j));
            endif;
        endfor;
        y0 += Y(i) * p;
    endfor;
    return;
endfunction

original_function(x0)

function y0=original_function(x0)
    y0 = exp(x0);
endfunction

delta_max(Y1, Y2)

function [delta, i] = delta_max(Y1,Y2)
    delta = 0;
    i = 1;
    n = size(Y1, 2);
    for j=1:n
        if abs(Y1(j) - Y2(j)) > delta
            delta = abs(Y1(j) - Y2(j));
            i = j;
        endif;
    endfor;
endfunction

I made some modifications with comments.

clear

function y0 = lagrange(X,Y,x0)
    y0 = 0;
    n = size(X, 2);
    for i = 1:n
        p = 1;
        for j = 1:n
            if (j != i)
##            Use .* for elementwise operating
              p = p.*(x0 - X(j)) / (X(i) - X(j));
            endif;
        endfor;
        y0 += Y(i) * p;
    endfor;
## I removed the "return" command
endfunction

function y0=original_function(x0)
    y0 = exp(x0);
endfunction

function [delta, i] = delta_max(Y1,Y2)
    delta = 0;
    i = 1;
    n = size(Y1, 2);
    for j=1:n
        if abs(Y1(j) - Y2(j)) > delta
            delta = abs(Y1(j) - Y2(j));
            i = j;
        endif;
    endfor;
endfunction

# You can call the functions directly by its name
X = [1, 2, 3]
Y = original_function(X)
h = (X(end) - X(1)) / (100 * (size(X, 2) - 1))
x_h = X(1):h:X(end)
y_h = original_function(x_h)
lagr = lagrange(X, Y, x_h)
[delta, x_delta_max] = delta_max(y_h, lagr);