Make Lagrange interpolation function faster?

Hi all,

I coded a Lagrange interpolation function for matrices, and calls it in my main program. only 4 points linear interpolation is considered. It runs for 1010400 times, which costs 397.3 seconds. Is there any ways to make it faster so that I could save some time running my main program?

Here is the function: Input coefficients from Lagrange polynomials, the coordinates (x, y) for the interpolated points, output lag_val gives the result of interpolation.

function [lag_val] = LagInterpolationOtptSingle(coeff, x, y)
  % when coeff are not consisted of square blocks but rectangular blocks,
  % goes to last case, where only 4 sample points linear case is considered.
  m = size(coeff, 1);
  n = size(coeff, 2);

    if n == 1

        if m == 3

            lag_val = coeff(1)*x+coeff(2)*y+coeff(3)*x^0;

        elseif m == 4

            lag_val = coeff(1)*x*y+coeff(2)*x+coeff(3)*y+coeff(4)*x^0;

        elseif m == 6

            lag_val = coeff(1)*x^2+coeff(2)*x*y+coeff(3)*y^2+...
                coeff(4)*x+coeff(5)*y+coeff(6)*x^0;

        elseif m == 8

            lag_val = coeff(1)*x^2*y+coeff(2)*x*y^2+...
                coeff(3)*x^2+coeff(4)*y^2+coeff(5)*x*y+coeff(6)*x+coeff(7)*y+...
                coeff(8)*x^0;

        elseif m == 9

            lag_val = coeff(1)*x^2*y^2+coeff(2)*x^2*y+coeff(3)*x*y^2+...
                coeff(4)*x^2+coeff(5)*y^2+coeff(6)*x*y+coeff(7)*x+coeff(8)*y+...
                coeff(9)*x^0;

        end

    elseif n>1

        ord = m/n;

        if ord == 3

            lag_val = coeff(1:n, :)*x+...
                coeff(2*n-n+1:2*n, :)*y+...
                coeff(3*n-n+1:3*n, :)*x^0;

        elseif ord == 4

            lag_val = coeff(1:n, :)*x*y+...
                coeff(2*n-n+1:2*n, :)*x+...
                coeff(3*n-n+1:3*n, :)*y+...
                coeff(4*n-n+1:4*n, :)*x^0;

        elseif ord == 6

            lag_val = coeff(1:n, 1:n)*x^2+...
                coeff(2*n-n+1:2*n, :)*x*y+...
                coeff(3*n-n+1:3*n, :)*y^2+...
                coeff(4*n-n+1:4*n, :)*x+...
                coeff(5*n-n+1:5*n, :)*y+...
                coeff(6*n-n+1:6*n, :)*x^0;

        elseif ord == 8

            lag_val = coeff(1:n, 1:n)*x^2*y+...
                coeff(2*n-n+1:2*n, :)*x*y^2+...
                coeff(3*n-n+1:3*n, :)*x^2+...
                coeff(4*n-n+1:4*n, :)*y^2+...
                coeff(5*n-n+1:5*n, :)*x*y+...
                coeff(6*n-n+1:6*n, :)*x+...
                coeff(7*n-n+1:7*n, :)*y+...
                coeff(8*n-n+1:8*n, :)*x^0;

        elseif ord == 9

            lag_val = coeff(1:n, 1:n)*x^2*y^2+...
                coeff(2*n-n+1:2*n, :)*x^2*y+...
                coeff(3*n-n+1:3*n, :)*x*y^2+...
                coeff(4*n-n+1:4*n, :)*x^2+...
                coeff(5*n-n+1:5*n, :)*y^2+...
                coeff(6*n-n+1:6*n, :)*x*y+...
                coeff(7*n-n+1:7*n, :)*x+...
                coeff(8*n-n+1:8*n, :)*y+...
                coeff(9*n-n+1:9*n, :)*x^0;

        else
            p = m/4;
            lag_val = coeff(1:p, :)*x*y+...
                coeff(2*p-p+1:2*p, :)*x+...
                coeff(3*p-p+1:3*p, :)*y+...
                coeff(4*p-p+1:4*p, :)*x^0; 
        end
    end