best way to evaluate a 2d function composed by a summation of functions

I need to plot this function:

where c is a vector of constants, and each Φ function is a polynomial expasion of the 4+i or 4+j order.

At the moment the functions are saved in two struct (dimensions 1xN and 1xM) such that in each dimention of one structure I have the function handle reletated to the x or the y coordinate.

I'm using the following lines to compute and to plot w :

%...
x = linspace(0,xmax,Nox);
y = linspace(0,ymax,Noy);
W = zeros(Nox,Noy);

for i = 1:N
    for j = 1:M
            W = W + c(i*j)*(trialx(i).PHI(x))'*(trialy(j).PHI(y));
    end
end

figure('Name','Displacement field','NumberTitle','off');
surf(x,y,W','EdgeColor','none'); daspect([100000 100000 1]);
xlabel('x [mm]','Fontsize',16); ylabel('y [mm]','Fontsize',16)
zlabel('y [mm]','Fontsize',16)

All the code is implemented with the symbolic toolbox and I'm struggling to plot w(x,y) with fsurf because I'm not able to write the correct summation of functions.

I need to raise the number of samples, N and M and it is quite slow;

Do you have any suggestions?