This is a MATLAB function for computing and plotting a cubic spline interpolation given a set of input data points. Additionally, "zoom_image" is another function that takes an image file path and a zoom factor as input, and outputs a zoomed version of the image, using cubic spline interpolation.

cubic_spline_interp.m

function [s] = cubic_spline_interp(X, Y, options)
    arguments
        X (1,:) {mustBeVector}
        Y (1,:) {mustBeVector}
        options.Plotting (1, 1) string = 'False'
    end

    % Calculation of differences between consecutive x-coordinates.
    n = length(X);
    h = zeros(n-1, 1);
    for i=1:n-1
        h(i) = X(i+1) - X(i);
    end
    
    % Construction of the coefficient matrix and the right-hand side vector
    % for solving the system of equations.
    A = zeros(n, n);
    b = zeros(n, 1);
    A(1,1) = 1;
    A(n,n) = 1;
    for i=2:n-1
        A(i, i-1) = h(i-1);
        A(i, i) = 2*(h(i-1) + h(i));
        A(i, i+1) = h(i);
        b(i) = 3 * ((Y(i+1) - Y(i))/h(i) - (Y(i) - Y(i-1))/h(i-1));
    end
    
    % Solution of the system of equations.
    c = A\b;

    % Calculation of coefficients of the first and second derivative for each
    % segment.
    b = zeros(n-1, 1);      % Coefficients of first derivative
    d = zeros(n-1, 1);      % Coefficients of second derivatives
    for i=1:n-1
        b(i) = (Y(i+1) - Y(i))/h(i) - h(i)*(c(i+1) + 2*c(i))/3;
        d(i) = (c(i+1) - c(i))/(3*h(i));
    end
    
    % Store the information about the spline
    s.X = X;
    s.Y = Y;
    s.h = h;
    s.c = c;
    s.d = d;
    s.b = b; 
    s.eval = @(x)eval_cubic_spline(x, s);

    % Plot the spline structure
    if isequal(options.Plotting, 'True')
        x_interp = linspace(min(X), max(X), 1000);

        % Evaluate the spline at each X in x_interp
        y_interp = s.eval(x_interp);
        
        % Plot the data points and the spline
        plot(x_interp, y_interp, 'b-', 'LineWidth', 1.5, 'Color', 'red');
        hold on;
        plot(X, Y, 'o', 'MarkerFaceColor', 'red', 'MarkerEdgeColor', 'black');
        hold off;
        legend('Cubic spline', 'Data Points');
        xlabel('X');
        ylabel('Y');
        grid on; box on; grid minor;
        set(gca,'FontSize',12,'LineWidth',1.5);
    end
end

function [y] = eval_cubic_spline(x, s)
    % Initialize a vector y to store the function values evaluated at x.
    num_x = length(x);
    y = zeros(num_x, 1);
    for i=1:num_x
        j = 1;
        while x(i) > s.X(j+1)
            j = j + 1;
        end
        
        x0 = s.X(j);
        y0 = s.Y(j);
        b = s.b(j);
        c = s.c(j);
        d = s.d(j);
        
        % The formula is used to evaluate the spline at point x.
        y(i) = y0 + b*(x(i)-x0) + c*(x(i)-x0)^2 + d*(x(i)-x0)^3;
    end
end

zoom_image.m

function zoom_image(image_path, z)
    % Read image
    I = imread(image_path);
    
    % Get image dimensions
    [H, W, ~] = size(I);
    
    % Generate x and y coordinates for original and zoomed image
    x = 1:W;
    y = 1:H;
    xi = linspace(1, W, z*W);
    yi = linspace(1, H, z*H);
    
    % Interpolate the image using cubic spline
    I = double(I);
    I_new = zeros(z*H, z*W, 3);
    for c = 1:3 % interpolate each color channel separately
        I_c = I(:, :, c);
        for i = 1:H
            % Evaluate the cubic spline at the new y-coordinate
            spline_y = cubic_spline_interp(y, I_c(i, :));
            I_new(i, :, c) = spline_y.eval(xi);
        end
        
        I_c_new = I_new(:, :, c);
        for j = 1:z*W
            % Evaluate the cubic spline at the new x-coordinate
            spline_x = cubic_spline_interp(x, I_c_new(:, j));
            I_new(:, j, c) = spline_x.eval(yi);
        end
    end
    
    % Convert image to uint8
    I_new = uint8(I_new);
    
    % Show image
    imshow(I_new);
    
    % Save image
    [path, name, ext] = fileparts(image_path);
    imwrite(I_new, fullfile(path, strcat(name, '_zoomed', ext)));
end