Error while solving for chebyshev

Hi. Im trying to run this code in MATLAB. but it gives me error. I dont know how to correct the error as I dont understand the error. Can anyone help me?

This is the error:

   Error using horzcat
Dimensions of matrices being concatenated are not consistent.

Error in Meanflowsolver (line 55)
M = sqrt(2/N).*[0.5.*ones(N+1,1),cos((J-1).*(K-1).*pi/(N)),0.5.*(-1).*(K-1)];
  end

    function [U, D2U] = Meanflowsolver(N);
% Function to transform the mean velocity profiles from collocation space
% to % Chebyshev space.
% Collocation values for U and U"
Ml = [0.000000000000000000; 0.000282959826495337; 0.001132694492131107;
0.002551775533010251; 0.010277021734211642; 0.023389962590386885;
0.042263021572698345; 0.067451726642569251; 0.099733544868301743;
0.140162787563543330; 0.190122654620963450; 0.251399474889032580;
0.326187775921186310; 0.416939393914891540; 0.525650219929556740;
0.651784704086493690; 0.787257935681726860; 0.908916569889709480;
0.981551413490585610; 0.999344391475489920; 0.999999792658503250;
0.999999999993113620; 1.000000000000001100; 0.004544508840776532;
0.014034420765414753; 0.029017609018465982; 0.049924151789873243;
0.077378438249061141; 0.112243586658369160; 0.155678082277931010;
0.209189502609298690; 0.274709727098769120; 0.354544401601672220;
0.451122324763552300; 0.565933832462586440; 0.696656823122624450;
0.831050446620130940; 0.940062426410669840; 0.992017489365257470;
0.999900693387071770; 0.999999998217581790; 0.999999999999882540;
1.000000000000000000; 0.007116964915322029; 0.018400836816274111;
0.035301822642907835; 0.058310865885997161; 0.088126108383279611;
0.125703085474088840; 0.172309413067983810; 0.229587783399974150;
0.299608078640860360; 0.384768565221385310; 0.487359576818445720;
0.608064048047086140; 0.742076148885554730; 0.872095613969664060;
0.964415132962072820; 0.997294784907640230; 0.999992311429745320;
1.000000000025288400; 1.000000000000046000; 1.000000000000000000;
1.000000000000000000];

M2 = [0.000000000000000000; -0.000000059285521856; -0.000000949996403597;
-0.000004821441696169; -0.000015291846444472; -0.000037503201616729;
-0.000078199640696901; -0.000145830968394501; -0.000250683486000673;
-0.000405040816390607; -0.000623378035768106; -0.000922593072477855;
-0.001322280034500118; -0.001845049868001658; -0.002516904508894791;
-0.003367671428042255; -0.004431506119054188; -0.005747470519947258;
-0.007360195409595532; -0.009320637882472041; -0.011686949375501972;
-0.014525221797035895; -0.017910587265962605; -0.021928356593549674;
-0.026675161512582653; -0.032260074630909921; -0.038805644088130142;
-0.046448035538733744; -0.055336399385716273; -0.065632659410805858;
-0.077508913993729250; -0.091141995976364895; -0.106702613879449790;
-0.124339165575672870; -0.144155443544363600; -0.166172178862819070;
-0.190273572482252460; -0.216133702460587970; -0.243112462357646470;
-0.270148350500162890; -0.295577792014319110; -0.317074012165833090;
-0.331629895227252390; -0.335578815616557000; -0.325204101062522450;
-0.297785052677842050; -0.253037686836590450; -0.194735124500390100;
-0.131333192176974560; -0.074194340142948109; -0.032961888571106240;
-0.010523218383097537; -0.002114337177937716; -0.000219088960191523;
-0.000008140072457660; -0.000000095887770295; 0.000000000944915165;
-0.000000001030438495; -0.000000000043537120; 0.000000000010812587;
0.000000000000473505; 0.000000000000014067; -0.000000000000000137;
0.000000000000000000];
% Use cosine transform to get Chebyshev values
Ml = [Ml; 1];
M2 = [M2; 0];
J = [2;N];
K = [1;N+1];
M = sqrt(2/N).*[0.5.*ones(N+1,1),cos((J-1).*(K-1).*pi/(N)),0.5.*(-1).*(K-1)];
U1 = sqrt(2/N)*M*Ml;
U2 = sqrt(2/N)*M*M2;
% Matrix multiplier for U and U"
for K = 1:N
for J = 1:N
J1 = abs(K-J)+l;
J2 = K+J;
U(K,J) = 0.5*U1(J1);
if (J2 <= N)
U(K,J) = U(K,J) + 0.5*U1(J2);
end
end
end
for K = 1 :N
for J = 1:N
J1 = abs(K-J)+l;
J2 = K+J;
D2U(K,J) = 0.5*U2(J1);
if (J2 <= N)
D2U(K,J) = D2U(K,J) + 0.5*U2(J2);
end
end
end
end