What's wrong in my code?

x = [1959:2019];
y = [315.98; 316.91; 317.64; 318.45; 318.99; 319.62; 320.04; 321.37; 322.18; 323.05; 324.62; 325.68; 326.32; 327.46; 329.68;
    330.18; 331.12; 332.04; 333.83; 335.40; 336.84; 338.75; 340.11; 341.45; 343.05; 344.66; 346.12; 347.43; 349.18; 351.57;
    353.12; 354.39; 355.61; 356.45; 357.10; 358.83; 360.82; 362.61; 363.73; 366.70; 368.38; 369.55; 371.14; 373.28; 375.80;
    377.52; 379.80; 381.90; 383.79; 385.59; 387.43; 389.90; 391.65; 393.86; 396.52; 398.64; 400.83; 404.22; 406.55; 408.52; 411.43];
    
%Perform least - squares regression
%Define basis function
phi = @(x) [1 exp(x) sin(x)];
%Creaate PHI matrix
PHI = cell2mat(arrayfun(phi, x, 'Uni', 0));
%Solve for W
W = (PHI'.*PHI)\(PHI'.*y);

%Display results
%Report resulting model
fprintf('y = %.8f + %.8f exp(x) + %.8f sin(x)\n', W);

%Finely spaced "x" values
xp = linspace(x(1), x(end), 100)';
PHI = cell2mat(arrayfun(phi, xp, 'Uni', 0));
%Predicted "y" at those "x"
yp = PHI*W;

%Calculate R2
yp = W(1) + W(2)*exp(x) + W(3)*sin(x);
R2 = 1 - sum((y - yp).^2)/sum((y - mean(y)).^2);
fprintf('R2: %.4f\n', R2);

[Walter Roberson]

You assign to yp twice in a row?

[dpb]

>> W = (PHI'.*PHI)\(PHI'.*y);
Matrix dimensions must agree. 
>> whos PHI x y
  Name      Size             Bytes  Class     Attributes

  PHI       1x183             1464  double              
  x          1x61              488  double              
  y         61x1               488  double              

>>

[dpb]

Also

>> phi(x(1))
ans =
    1.0000       Inf   -0.9766
>> 

exp(x) --> Inf for all values of x