error using symengine is generated while getting the derivative of equation of 8 degree with variable constants

I am trying to get the slope of equation (y3c=p1*xs^8 + p2*xs^7 + p3*xs^6 + p4*xs^5 + p5*xs^4 + p6*xs^3 + p7*xs^2 + p8*xs + p9),

where,xs=(x3s-mean)./dev; and xs x3s are symbols.

The problem is that the constant p1 to p9 ,mean and dev are arrays.

how should i modify the code to work without the following error

Error using symengine

Matrix must be square.

Error in sym/privBinaryOp (line 1034)

Csym = mupadmex(op,args{1}.s, args{2}.s, varargin{:});

Error in ^ (line 339)

B = privBinaryOp(A, p, 'symobj::mpower');

Error in trailing_edge_2 (line 79)

y3c=p1*xs^8 + p2*xs^7 + p3*xs^6 + p4*xs^5 + p5*xs^4 + p6*xs^3 + p7*xs^2 + p8*xs + p9;

clc
clear
Excel= xlsread('cf.xlsx','cf') ;
x=Excel(:,1);%x
y=Excel(:,2);%y

%%%%%%%%%%%%trailing edge(551:614,1:36))%%%%%%%%%%%
x3=[x(551:614);x(1:36)];
y3=[y(551:614);y(1:36)];
x3=x3';
y3=y3';
for i=1:length(x3)
    if i<=33

        mean(i)=-0.0005415;dev(i)=0.0007397;
       p1(i) =  -8.172e-05;p2(i) =  -0.0004409 ; p3(i) =  -0.0006484;  
       p4(i) =    0.000167;p5(i) =   0.0008287;  p6(i) =    1.85e-05;
       p7(i) =  -0.0006773;p8(i) =   -0.001281;  p9(i) =    0.004286;  
    elseif (i>=34) && (i<=39)
 
       mean(i)=7.32e-05;dev(i)=0.0001512;
       p1(i) =   5.645e-05  ;p2(i) =    0.000194  ;p3(i) =   7.821e-05  ;
       p4(i) =  -0.0002781  ;p5(i) =   -0.000149  ;p6(i) =   0.0002006  ;
       p7(i) =   0.0002051  ;p8(i)=   0.0005569   ;p9(i) =   0.0002695  ;          
    else

        mean(i)=-0.002376; dev(i)=0.002098;
       p1(i) =   1.409e-05  ;p2(i) =   8.773e-05  ;p3(i) =   0.0001679  ;
       p4(i) =   5.782e-05  ;p5(i) =  -8.276e-05  ;p6(i) =    0.000105  ;
       p7(i) =   0.0006627  ;p8(i) =    0.001226  ;p9(i) =   -0.002926  ;
    end

    x3p(i)=(x3(i)-mean(i))/dev(i);
    y3p(i)=p1(i)*x3p(i)^8 + p2(i)*x3p(i)^7 + p3(i)*x3p(i)^6 + p4(i)*x3p(i)^5 + p5(i)*x3p(i)^4 + p6(i)*x3p(i)^3 + p7(i)*x3p(i)^2 + p8(i)*x3p(i) + p9(i);%curve plotted using equ.

end
diff22=max(y3-y3p);
syms xs x3s
xs=(x3s-mean)./dev;
y3c=p1*xs^8 + p2*xs^7 + p3*xs^6 + p4*xs^5 + p5*xs^4 + p6*xs^3 + p7*xs^2 + p8*xs + p9;
dys=diff(y3c);
slope=double(vpa(subs(dys,x3s,x3)));
theta=atand(slope); %%% slope array