Save outputs of nested for loops in MATLAB

I have the following codes which I wish to have an output matrix Rpp of 10201 rows and 3 columns. I run this code (which takes a bit long) then I check the matrix size of Rpp and I see a vector of 1 row and 3 columns, I tried so many things I couldn't find any proper way. The logic of the codes is to take the 6 values (contain 4 constant values and 2 variable values (chosen from 101 values)) and make the calculation for 3 different i1 and store every output vector of 3 in a matrix with (101*101 (pairs of those 2 variable values)) rows and 3 (for each i1) columns.

I appreciate your help

      Vp1=linspace(3000,3500,101);
      Vp2=3850;
      rho1=2390;
      rho2=2510;
      Vs1=linspace(1250,1750,101);
      Vs2=2000;
      i1=[10 25 40];

      Rpp = zeros(length(Vp1)*length(Vs1),length (i1));

      for n=1:length(Vp1)*length(Vs1)
      for m=1:length (i1)
      for l=1:length(Vp1)
          for k=1:length(Vs1)

             p=sin(i1)/Vp1(l);
             i2=asin(p*Vp2);
             j1=asin(p*Vs1(k));
             j2=asin(p*Vs2);
             a=rho2*(1-2*Vs2^2*p.^2)-rho1*(1-2*Vs1(k).^2*p.^2);
             b=rho2*(1-2*Vs2^2*p.^2)+2*rho1*Vs1(k)^2*p.^2;
             c=rho1*(1-2*Vs1(k)^2*p.^2)+2*rho2*Vs2^2*p.^2;
             d=2*(rho2*Vs2^2-rho1*Vs1(k)^2);
             E=b.*cos(i1)./Vp1(l)+c.*cos(i2)/Vp2;
             F=b.*cos(j1)./Vs1(k)+c.*cos(j2)/Vs2;
             G=a-d*(cos(i1)/Vp1(l)).*(cos(j2)/Vs2);
             H=a-d*(cos(i2)/Vp2).*(cos(j1)/Vs1(k));
             D=E.*F+G.*H.*p.^2;
             Rpp=((b.*(cos(i1)/Vp1(l))-c.*cos((i2)/Vp2)).*F-(a+d*((cos(i1)/Vp1(l))).*(cos(j2)/Vs2)).*H.*p.^2)./D
          end
      end
      end
      end

Answer 1

[Walter Roberson]

      Rpp = zeros(length(Vp1), length(Vs1), length(i1));

      for l=1:length(Vp1)
          for k=1:length(Vs1)
             p=sin(i1)/Vp1(l);
             i2=asin(p*Vp2);
             j1=asin(p*Vs1(k));
             j2=asin(p*Vs2);
             a=rho2*(1-2*Vs2^2*p.^2)-rho1*(1-2*Vs1(k).^2*p.^2);
             b=rho2*(1-2*Vs2^2*p.^2)+2*rho1*Vs1(k)^2*p.^2;
             c=rho1*(1-2*Vs1(k)^2*p.^2)+2*rho2*Vs2^2*p.^2;
             d=2*(rho2*Vs2^2-rho1*Vs1(k)^2);
             E=b.*cos(i1)./Vp1(l)+c.*cos(i2)/Vp2;
             F=b.*cos(j1)./Vs1(k)+c.*cos(j2)/Vs2;
             G=a-d*(cos(i1)/Vp1(l)).*(cos(j2)/Vs2);
             H=a-d*(cos(i2)/Vp2).*(cos(j1)/Vs1(k));
             D=E.*F+G.*H.*p.^2;
             Rpp(l, k, :) = ((b.*(cos(i1)/Vp1(l))-c.*cos((i2)/Vp2)).*F-(a+d*((cos(i1)/Vp1(l))).*(cos(j2)/Vs2)).*H.*p.^2)./D;
          end
      end

Note: it would be more efficient if you looped in the other order:

      Rpp = zeros(length(i1), length(Vs1), length(Vp)1);

      for k=1:length(Vs1)
          for l=1:length(Vp1)
             p=sin(i1)/Vp1(l);
             i2=asin(p*Vp2);
             j1=asin(p*Vs1(k));
             j2=asin(p*Vs2);
             a=rho2*(1-2*Vs2^2*p.^2)-rho1*(1-2*Vs1(k).^2*p.^2);
             b=rho2*(1-2*Vs2^2*p.^2)+2*rho1*Vs1(k)^2*p.^2;
             c=rho1*(1-2*Vs1(k)^2*p.^2)+2*rho2*Vs2^2*p.^2;
             d=2*(rho2*Vs2^2-rho1*Vs1(k)^2);
             E=b.*cos(i1)./Vp1(l)+c.*cos(i2)/Vp2;
             F=b.*cos(j1)./Vs1(k)+c.*cos(j2)/Vs2;
             G=a-d*(cos(i1)/Vp1(l)).*(cos(j2)/Vs2);
             H=a-d*(cos(i2)/Vp2).*(cos(j1)/Vs1(k));
             D=E.*F+G.*H.*p.^2;
             Rpp(:, l, k) = ((b.*(cos(i1)/Vp1(l))-c.*cos((i2)/Vp2)).*F-(a+d*((cos(i1)/Vp1(l))).*(cos(j2)/Vs2)).*H.*p.^2)./D;
          end
      end