Matlab - Coeficiente Correlacion r2

Tengo el siguiente programa en MATLAB7

clc;clear;
datos26; %Archivo de ingreso de datos
P1=0.000;P2=0.000;P3=0.000;P4=0.000;P5=0.000;P6=0.000;
n=length(H2G);
Co=0;
Rx=0;
opcion=0;sg=1;ajuste=1;xn=H42;yn=H2G; %Condiciones iniciales
% CALCULAMOS H26
%H26=H5C2*RES2^P1/RES4^P2/H39^P3*H40^P4/H45^P5*H47^P6
H26=H5C2.*RES2.^P1./RES4.^P2./H39.^P3.*H40.^P4./H45.^P5.*H47.^P6;
% CALCULAMOS H36 -> H36=H26
H36=H26;
% CALCULAMOS H42 -> H42 = H36/H3C*100
H42=H36./H3C.*100;
%Grado del Polinomio
m=input(' Grado del polinomio = ');
if m>n
ajuste=0;
end

if ajuste==0
fprintf('\n No se puede realizar el ajuste \n');
else
for i=1:m+1
for j=1:m+1
sx(i,j)=sum(xn.^(i+j-2));
end
sy(i)=sum(yn.*xn.^(i-1));
end

c=sx\sy';
xx=linspace(min(H42),max(H42));
cn=flipud(c);
ya=polyval(cn,H42);
yy=polyval(cn,xx);
%Calculo coeficiente correlacion
st=sum((H2G-mean(H2G)).^2);
sr=sum((H2G-ya).^2);
r=sqrt((st-sr)/st);
% Hacemos el Bucle de While y For.
while r <= 0.80
for P6=-15.00:+0.05:15.00
for P5=-15.00:+0.05:15.00
for P4=-15.00:+0.05:15.00
for P3=-15.00:+0.05:15.00
for P2=-15.00:+0.05:15.00
for P1=-15.00:+0.05:15.00

% CALCULAMOS H26
% H26=H5C2*RES2^P1/RES4^P2/H39^P3*H40^P4/H45^P5*H47^P6
H26=H5C2.*RES2.^P1./RES4.^P2./H39.^P3.*H40.^P4./H45.^P5.*H47.^P6;
% CALCULAMOS H36 -> H36=H26
H36=H26;
% CALCULAMOS H42 -> H42 = H36/H3C*100
H42=H36./H3C.*100;

for i=1:m+1
for j=1:m+1
sx(i,j)=sum(xn.^(i+j-2));
end
sy(i)=sum(yn.*xn.^(i-1));
end

%Polyval
c=sx\sy';
xx=linspace(min(H42),max(H42));
cn=flipud(c);
ya=polyval(cn,H42);
yy=polyval(cn,xx);
%Calculo coeficiente correlacion
st=sum((H2G-mean(H2G)).^2);
sr=sum((H2G-ya).^2);
r=sqrt((st-sr)/st);

if r >= Rx
Rx=r;
%Mostramos resultados
fprintf('\n r= %g ',Rx);
fprintf('P1= %g ',P1);
fprintf('P2= %g ',P2);
fprintf('P3= %g ',P3);
fprintf('P4= %g ',P4);
fprintf('P5= %g ',P5);
fprintf('P6= %g ',P6);
fprintf('Co= %g ',Co);
end
%Incrementamos el numero de combinaciones
Co=Co+1;
end
end
end
end
end
end
end

fprintf('\n Coeficiente de correlacion: r = %g \n',r);
fprintf('\n Presione cualquier tecla para ver la grafica del ajuste \n');
pause
plot(xx,yy,H42,H2G,'.r');
pause
close all
end

El fichero datos es de este estilo...
H42=[38.90035 38.85169 38.82986 38.79698 38.67423 40.11880 40.13742 40.15181 40.18371 40.23192 39.85407 39.95980 39.94731];
H13=[36.00694 36.00694 36.00694 36.00694 36.00694 37.27668 37.27668 37.27668 37.27668 37.27668 37.49313 37.49313 37.49313 ];
H2G=[73 73 73 73 73 84 84 84 84 84 86 86 86 93 93 93 93 93];
H5C2=[24.29895 24.60516 24.73010 24.91754 25.44995 17.72882 17.79633 17.86481 17.84717 17.77112 22.52222 22.34137 23.01654];
RES2=[1.02483 1.02562 1.02599 1.02656 1.02895 1.00841 1.00778 1.00724 1.00653 1.00560 1.00408 1.00185 1.00098];
RES4=[1.09462 1.09034 1.08622 1.08350 1.08008 1.20282 1.20000 1.19540 1.19311 1.19236 1.04332 1.03562 1.01972];
H39=[1.04261 1.03431 1.03116 1.02652 1.01644 1.26672 1.25945 1.25272 1.24974 1.24868 1.02488 1.02200 0.99330];
H40=[1.01669 1.01756 1.01798 1.01865 1.02149 1.01232 1.01140 1.01059 1.00990 1.00919 0.99616 0.99405 0.99252];
H45=[1.00400 1.00416 1.00457 1.00480 1.00592 1.00672 1.00602 1.00566 1.00532 1.00488 0.99806 0.99784 0.99591];
H47=[1.01747 1.01844 1.01891 1.01965 1.02271 1.01246 1.01163 1.01090 1.01028 1.00966 0.99748 0.99550 0.99404];
H3C=[62.46460 63.33099 63.68835 64.22546 65.80597 44.19080 44.33850 44.49316 44.41394 44.17169 56.51172 55.90961 57.61725];

Cada una de las matrices son de 400 datos.
Intento ver el r (coef.Corr.) mejor , modificando los datos de P1,P2,P3,P4,P5 y P6.
Todo funciona si m=2, es decir si el polinomio es de Grado 2. Si pongo m=3 (de Grado 3), el resultado obtenido no coincide con Excel o StartGraphics.
Que tengo que modificar para conseguir r2, para Grado 3.
Gracias. Y Feliz Año 2010.