Need help to overcome this error 'Dimensions of matrices being concatenated are not consistent.' Any idea/comment is appreciated!

I have the following code:

      phi1_12 = 0;
      phi1_13 = 0;

      A = pi;
      a = cos(A/2);
      b = -1i*sin(A/2);
      ai = cos(A/2);
      bi = -b;

      S = @(theta)sin(theta);
      C = @(theta)cos(theta);

      zeta = 1;

      ep1_12 = exp(1i*phi1_12);
      ep1_13 = exp(1i*phi1_13);

      ep1_12c = exp(-1i*phi1_12);
      ep1_13c = exp(-1i*phi1_13);

       ep2_12 = @(phi2_12)exp(1i*phi2_12);
       ep2_13 = @(phi2_13)exp(1i*phi2_13);
      % 
       ep2_12c = @(phi2_12)exp(-1i*phi2_12);
       ep2_13c = @(phi2_13)exp(-1i*phi2_13);
      % 
       ep3_12 = @(phi3_12)exp(1i*phi3_12);
       ep3_13 = @(phi3_13)exp(1i*phi3_13);
      % 
       ep3_12c = @(phi3_12)exp(-1i*phi3_12);
       ep3_13c = @(phi3_13)exp(-1i*phi3_13);
      % 
      ep4_12 = @(phi4_12)exp(1i*phi4_12);
      ep4_13 = @(phi4_13)exp(1i*phi4_13);
      % 
       ep4_12c = @(phi4_12)exp(-1i*phi4_12);
       ep4_13c = @(phi4_13)exp(-1i*phi4_13);
      % 
       ep5_12 = @(phi5_12)exp(1i*phi5_12);
       ep5_13 = @(phi5_13)exp(1i*phi5_13);
      % 
       ep5_12c = @(phi5_12)exp(-1i*phi5_12);
       ep5_13c = @(phi5_13)exp(-1i*phi5_13);

       U1 = @(theta)[a,b.*ep1_12.*C(theta),b.*ep1_13.*S(theta);-bi.*ep1_12c.*C(theta),ai.*C(theta).^2 + zeta.*S(theta).^2,(ai - zeta).*exp(-1i.*(phi1_12-phi1_13)).*S(theta).*C(theta);-bi.*ep1_13c.*S(theta),(ai - zeta).*exp(1i.*(phi1_12-phi1_13)).*S(theta).*C(theta),ai.*S(theta).^2 + zeta.*C(theta).^2];
       U2 = @(theta,phi2_12,phi2_13)[a,b.*ep2_12(phi2_12).*C(theta),b.*ep2_13(phi2_13).*S(theta);-bi.*ep2_12c(phi2_12).*C(theta),ai.*C(theta).^2 + zeta.*S(theta).^2,(ai - zeta).*exp(-1i.*(phi2_12-phi2_13)).*S(theta).*C(theta);-bi.*ep2_13c(phi2_13).*S(theta),(ai - zeta).*exp(1i.*(phi2_12-phi2_13)).*S(theta).*C(theta),ai.*S(theta).^2 + zeta.*C(theta).^2];
       U3 = @(theta,phi3_12,phi3_13)[a,b.*ep3_12(phi3_12).*C(theta),b.*ep3_13(phi3_13).*S(theta);-bi.*ep3_12c(phi3_12).*C(theta),ai.*C(theta).^2 + zeta.*S(theta).^2,(ai - zeta).*exp(-1i.*(phi3_12-phi3_13)).*S(theta).*C(theta);-bi.*ep3_13c(phi3_13).*S(theta),(ai - zeta).*exp(1i.*(phi3_12-phi3_13)).*S(theta).*C(theta),ai.*S(theta).^2 + zeta.*C(theta).^2];
       U4 = @(theta,phi4_12,phi4_13)[a,b.*ep4_12(phi4_12).*C(theta),b.*ep4_13(phi4_13).*S(theta);-bi.*ep4_12c(phi4_12).*C(theta),ai.*C(theta).^2 + zeta.*S(theta).^2,(ai - zeta).*exp(-1i.*(phi4_12-phi4_13)).*S(theta).*C(theta);-bi.*ep4_13c(phi4_13).*S(theta),(ai - zeta).*exp(1i.*(phi4_12-phi4_13)).*S(theta).*C(theta),ai.*S(theta).^2 + zeta.*C(theta).^2];
       U5 = @(theta,phi5_12,phi5_13)[a,b.*ep5_12(phi5_12).*C(theta),b.*ep5_13(phi5_13).*S(theta);-bi.*ep5_12c(phi5_12).*C(theta),ai.*C(theta).^2 + zeta.*S(theta).^2,(ai - zeta).*exp(-1i.*(phi5_12-phi5_13)).*S(theta).*C(theta);-bi.*ep5_13c(phi5_13).*S(theta),(ai - zeta).*exp(1i.*(phi5_12-phi5_13)).*S(theta).*C(theta),ai.*S(theta).^2 + zeta.*C(theta).^2];
      % 
      % 
      % 
      % 
       U = @(theta,phi2_13,phi2_12,phi3_12,phi3_13,phi4_12,phi4_13,phi5_12,phi5_13)U5(theta,phi5_12,phi5_13)*U4(theta,phi4_12,phi4_13)*U3(theta,phi3_12,phi3_13)*U2(theta,phi2_12,phi2_13)*U1(theta);

       sel = @(U,r,c)U(r,c);  % indexing the U(2,1) matrix element
       U21 = @(theta,phi2_13,phi2_12,phi3_12,phi3_13,phi4_12,phi4_13,phi5_12,phi5_13)sel(U(theta,phi2_13,phi2_12,phi3_12,phi3_13,phi4_12,phi4_13,phi5_12,phi5_13),2,1);

      N = 5;
      t = pi/4;
      % 
      % x = zeros(20);

      for i = 1:N
          phi1 = i*pi/N;
          for j = 1:N
              phi2 = j*pi/N;
               for k = 1:N
                   phi3 = k*pi/N;
                   for l = 1:N
                       phi4 = l*pi/N;
                        for m = 1:N
                            phi5 = m*pi/N;
                            for n = 1:N
                                phi6 = n*pi/N;

                                for o = 1:N
                                    phi7 = o*pi/N;
                                    for p = 1:N
                                        phi8 = p*pi/N;

                                        J(i,j,k,l,m,n,o,p) = abs((1/t)*integral(@(theta)real(U21(theta,phi1,phi2,phi3,phi4,phi5,phi6,phi7,phi8)),0,t) - real(U21(0,phi1,phi2,phi3,phi4,phi5,phi6,phi7,phi8))) + abs((1/t)*integral(@(theta)abs(U21(theta,phi1,phi2,phi3,phi4,phi5,phi6,phi7,phi8)),0,t) - 1);

                                    end 

                                end 

                            end 

                        end

                   end

               end

          end

      end

I already added '.*' and '.^' instead of * and ^ but still I am confused why the matrices are not consistent. I had checked them individually like for U1,U2... they do work!.. but there's problem while integrating it inside the for loop.

Answer 1

[Walter Roberson]

http://www.mathworks.com/help/matlab/ref/integral.html

For scalar-valued problems, the function y = fun(x) must accept a vector argument, x, and return a vector result, y

Your U1 through U5 functions are constructing arrays of results which do not take into account that theta will not be a scalar.

[JayashP]

Thank you for the answer. Can you please elaborate a bit on how to overcome this?

[Walter Roberson]

       U = @(th, phi2_13, phi2_12, phi3_12, phi3_13, phi4_12, phi4_13, phi5_12, phi5_13) arrayfun( @(theta) U5(theta,phi5_12,phi5_13) * U4(theta,phi4_12,phi4_13) * U3(theta,phi3_12,phi3_13) * U2(theta,phi2_12,phi2_13) * U1(theta), th) ;