Please, I would appreciate any assistance on debugging the attached code.

The first section of the code is okay, but I am having it difficult solving for uu and ud in the attached code. I want to continue solving for the values of uu and ud as long the conditions in the while loops are valid. Pls, see the attached file.

This is the pseudo code:

for j=1:N
    k=1;
    uu(k,j)=0;
    Lou(k,j)=2;
    p(k,j)=6000;
    pQ(k,j)=asin(1+Lou(k,j));
    Lh2(j)=200;
    k22(k,j)=Lh2(j)./Lou(k,j);

      if p(k,j)>pQ(k,j)
          while Lou(k,j)<Lh2
              k22f(k,j)=Lh2(j)./(Lou(k,j)+uus);
              funLu(k,j)=asin(1+k22f(k,j));
              testfunLu=@(x)eval(subs(funLu(k,j),uus,x))-p(k,j);
              uuf(k,j)=fsolve(testfunLu,55);
              uu(k+1,j)=uuf(k,j);
              Lou(k+1,j)=Lou(k,j)+uu(k+1,j);
              k22(k+1,j)=Lh2(j)./(Lou(k+1,j)+uu(k+1,j));
              p(k+1,j)=p(k,j);
              k=k+1
          end
      end
  end

[Chad Greene]

Nothing attached.

[Isaac]

This is the file. Thanks

[Isaac]

This is the part of the code, that has bug:

if p(k,j)>=pQu1(k,j)| p(k,j)>=pQd1(k,j) %displacements within the host formation (Case 1)

                      while Lou(k,j)<Lh2 & Lod(k,j)<Lh2
                          %Numerically Solve for ud at t=k+1 using Newton Raphson
                      funLu(k,j)=KICfm(j)./sqrt(pi*(Lou(k,j)+uus))+S2+0.5*grav*(Lou(k,j)+uus)*rhom; %up-tip
                      uuf{k,j}=real(double(solve(funLu(k,j)==p(k,j),uus)));
                      uu(k+1,j)=uuf{k,j}(2);
                          if uu(k+1,j)<0
                              uu(k+1,j)=0;
                          end
                     funLd(k,j)=KICfm(j)./sqrt(pi*(Lod(k,j)+uds))+S2-0.5*grav*(Lod(k,j)+uds)*rhom; %down-tip
                     udf{k,j}=real(double(solve(funLd(k,j)==p(k,j),uds)));
                     ud(k+1,j)=udf{k,j}(1);
                          if ud(k+1,j)<0
                              ud(k+1,j)=0;
                          end
                     Lou(k+1,j)=Lou(k,j)+uu(k+1,j);
                     pQu1(k+1,j)=KICfm(j)./sqrt(pi*Lou(k+1,j))+S2+0.5*grav*Lou(k+1,j)*rhom;
                     Lod(k+1,j)=Lod(k,j)+ud(k+1,j);
                     pQd1(k+1,j)=KICfm(j)./sqrt(pi*Lod(k+1,j))+S2-0.5*grav*Lod(k+1,j)*rhom;

                     p(k+1,j)=p(k,j);

                      %Case 2
                      pQd2(k+1,j)= pQd2(k,j);
                      pQu3(k+1,j)=pQu3(k,j);
                      k=k+1;

                      end

              if p(k,j)>pQd2(k,j)& Lod(k,j)>Lh2 & Lou(k,j)<Lh2 %case-2: up-tip within host formation and down-tip outside
                      Lod(k,j)=Lh2(1);
                      %Lou(k,j)=Lh2(1);

                      while (Lou(k,j)<=Lh2 & Lh2<Lod(k,j)<Lh3 )
                      %Calculating the displacements
                      %lower tip
                      %down-tip critical pressure function at Q
                      k22f(k,j)=Lh2(j)./(Lod(k,j)+uds);
                      funLd(k,j)=(sqrt(pi)*KICd(j)./(Lod(k,j)+uds)-grav*(Lod(k,j)+uds)*pi*rhom-sqrt((1+k22f(k,j))/(1-k22f(k,j)))*S2+k22f(k,j).*sqrt((1+k22f(k,j))/(1-k22f(k,j)))*S2+sqrt(1-k22f(k,j).^2)*S3+pi*S3+grav*(Lod(k,j)+uds)*rhom*asin(sqrt((1-k22f(k,j))/2))+(grav*(Lod(k,j)+uds)*rhom+2*S2-2*S3).*asin(sqrt((1+k22f(k,j))/2)))./(-sqrt((1+k22f(k,j))/(1-k22f(k,j)))+k22f(k,j).*sqrt((1+k22f(k,j))/(1-k22f(k,j)))+sqrt(1-k22f(k,j).^2)+pi);%down-tip
                      udf{k,j}=real(double(solve(funLd(k,j)==p(k,j),uds)));
                      ud(k+1,j)=udf{k,j}(2);
                      if ud(k+1,j)<0
                          ud(k+1,j)=0;
                      end            

                      %upper tip
                      k22(k,j)=Lh2(j)./(Lo(k,j)+ud(k,j));
                      funLu(k,j)=-(2*KICfm(j)*sqrt(pi)+sqrt((Lou(k,j)+uus)).*(grav*(Lou(k,j)+uus)*rhom*pi+2*(sqrt(1-k22(k,j).^2)*S2-(sqrt((1-k22(k,j))/(1+k22(k,j)))+k22(k,j).*sqrt((1-k22(k,j))/(1+k22(k,j)))-pi)*S3))+4*sqrt((Lou(k,j)+uus))*(S2-S3).*asin(sqrt((1+k22(k,j))/2)))./(2*sqrt((Lou(k,j)+uus)).*(sqrt((1-k22(k,j))/(1+k22(k,j)))+k22(k,j).*sqrt((1-k22(k,j))/(1+k22(k,j)))-sqrt(1-k22(k,j).^2)-pi)); %up-tip
                      uuf{k,j}=real(double(solve(funLu(k,j)==p(k,j),uus)));
                      uu(k+1,j)=uuf{k,j}(2);
                          if uu(k+1,j)<0
                              uu(k+1,j)=0;
                          end
                      Lod(k+1,j)=Lod(k,j)+ud(k+1,j);
                      k22(k+1,j)=Lh2(j)./(Lod(k+1,j)+ud(k+1,j));
                      pQd2(k+1,j)=(sqrt(pi)*KICd(j)./Lod(k+1,j)-grav*Lod(k+1,j)*pi*rhom-sqrt((1+k22(k+1,j))/(1-k22(k+1,j)))*S2+k22(k+1,j).*sqrt((1+k22(k+1,j))/(1-k22(k+1,j)))*S2+sqrt(1-k22(k+1,j).^2)*S3+pi*S3+grav*Lod(k+1,j)*rhom*asin(sqrt((1-k22(k+1,j))/2))+(grav*Lod(k+1,j)*rhom+2*S2-2*S3).*asin(sqrt((1+k22(k+1,j))/2)))./(-sqrt((1+k22(k+1,j))/(1-k22(k+1,j)))+k22(k+1,j).*sqrt((1+k22(k+1,j))/(1-k22(k+1,j)))+sqrt(1-k22(k+1,j).^2)+pi);%down-tip
                      Lou(k+1,j)=Lou(k,j)+uu(k+1,j);
                      pQu2(k+1,j)=-(2*KICfm(j)*sqrt(pi)+sqrt(Lou(k+1,j)).*(grav*Lou(k+1,j)*rhom*pi+2*(sqrt(1-k22(k+1,j).^2)*S2-(sqrt((1-k22(k+1,j))/(1+k22(k+1,j)))+k22(k+1,j).*sqrt((1-k22(k+1,j))/(1+k22(k+1,j)))-pi)*S3))+4*sqrt(Lou(k+1,j))*(S2-S3).*asin(sqrt((1+k22(k+1,j))/2)))./(2*sqrt(Lou(k+1,j)).*(sqrt((1-k22(k+1,j))/(1+k22(k+1,j)))+k22(k+1,j).*sqrt((1-k22(k+1,j))/(1+k22(k+1,j)))-sqrt(1-k22(k+1,j).^2)-pi)); %up-tip

                      p(k+1,j)=p(k,j);

                      %Case 3
                      pQu3(k+1,j)= pQu3(k,j);

                      k=k+1;

                      end
               elseif p(k,j)>pQu3(k,j)& Lou(k,j)>Lh2 & Lod(k,j)<Lh2 %case-3: down-tip within host formation and up-tip outside      

                      %Lod(k,j)=Lh2(1);
                      Lou(k,j)=Lh2(1);
                       while(Lod(k,j)<=Lh2 & Lh2<Lou(k,j)<Lh3) 
                     %Calculating the displacements

                      %upper tip
                      k13f(k,j)=Lh2(j)./(Lou(k,j)+uus);
                      funLu(k,j)=(2*KICu(j)*sqrt(pi)+sqrt((Lou(k,j)+uus)).*(grav*(Lou(k,j)+uus)*rhom*pi+2*(sqrt(1-k13f(k,j).^2)*S1+(-sqrt((1+k13f(k,j))./(1-k13f(k,j)))+k13f(k,j).*sqrt((1+k13f(k,j))./(1-k13f(k,j)))+pi)*S2))+4*sqrt((Lou(k,j)+uus))*(S1-S2).*asin(sqrt((1-k13f(k,j))/2)))./(2*sqrt((Lou(k,j)+uus)).*(-sqrt((1+k13f(k,j))./(1-k13f(k,j)))+k13f(k,j).*sqrt((1+k13f(k,j))./(1-k13f(k,j)))+sqrt(1-k13f(k,j).^2)+pi)); %up-tip
                      testfunLu=@(x)eval(subs(funLu(k,j),uus,x))-p(k,j);
                      %uuf{k,j}=real(double(fsolve(funLu(k,j)==p(k,j),uus)));
                      uuf(k,j)=fsolve(testfunLu,55);
                      uu(k+1,j)=uuf(k,f);
                      if uu(k+1,j)<0
                          uu(k+1,j)=0;
                      end
                      %lower tip
                      k13(k,j)=Lh2(j)./(Lou(k,j)+uu(k,j));
                      funLd(k,j)=0.5./(sqrt(Lod(k,j))*pi).*(2*KICfm(j)*sqrt(pi)+sqrt((Lod(k,j)+uds)).*(-grav*(Lod(k,j)+uds)*pi*rhom+2*(-sqrt((1-k13(k,j))/(1+k13(k,j))).*(1+k13(k,j))*S1+(sqrt((1-k13(k,j))/(1+k13(k,j)))+k13(k,j).*sqrt((1-k13(k,j))/(1+k13(k,j)))+pi)*S2))+4*sqrt((Lod(k,j)+uds))*(S1-S2).*asin(sqrt((1-k13(k,j))/2))); %down-tip
                      testfunLd=@(x)eval(subs(funLd(k,j),uds,x))-p(k,j);
                      udf(k,j)=fsolve(testfunLd,55);
                      ud(k+1,j)=udf(k,f);
                      if ud(k+1,j)<0
                          ud(k+1,j)=0;
                      end

                      Lou(k+1,j)=Lou(k,j)+uu(k+1,j);
                      k13(k+1,j)=Lh2(j)./(Lou(k+1,j)+uu(k+1,j));
                      pQu3(k+1,j)=(2*KICu(j)*sqrt(pi)+sqrt((Lou(k+1,j))).*(grav*(Lou(k+1,j))*rhom*pi+2*(sqrt(1-k13(k+1,j).^2)*S1+(-sqrt((1+k13(k+1,j))/(1-k13(k+1,j)))+k13(k+1,j).*sqrt((1+k13(k+1,j))/(1-k13(k+1,j)))+pi)*S2))+4*sqrt((Lou(k+1,j)))*(S1-S2).*asin(sqrt((1-k13(k+1,j))/2)))./(2*sqrt((Lou(k+1,j))).*(-sqrt((1+k13(k+1,j))/(1-k13(k+1,j)))+k13(k+1,j).*sqrt((1+k13(k+1,j))/(1-k13(k+1,j)))+sqrt(1-k13(k+1,j).^2)+pi)); %up-tip
                      Lod(k+1,j)=Lod(k,j)+ud(k+1,j);
                      pQd3(k+1,j)=0.5./(sqrt(Lod(k+1,j))*pi).*(2*KICfm(j)*sqrt(pi)+sqrt(Lod(k+1,j)).*(-grav*Lod(k+1,j)*pi*rhom+2*(-sqrt((1-k13(k+1,j))/(1+k13(k+1,j))).*(1+k13(k+1,j))*S1+(sqrt((1-k13(k+1,j))/(1+k13(k+1,j)))+k13(k+1,j).*sqrt((1-k13(k+1,j))/(1+k13(k+1,j)))+pi)*S2))+4*sqrt(Lod(k+1,j))*(S1-S2).*asin(sqrt((1-k13(k+1,j))/2))); %down-tip
                       p(k+1,j)=p(k,j);

                      %Case 4
                      pQu4(k+1,j)= pQu4(k,j);
                      pQd4(k+1,j)= pQd4(k,j);

                      k=k+1;
                       end
              end

      end

end

[Chad Greene]

Yikes. I can tell right away that you're unlikely to get a response because there's too much to sift through. Please trim your code down to a MWE and clearly describe what you mean when you say you're having difficulty.

[Isaac]

This is a summary:

for j=1:N
  k=1
 uu(k,j)=0 %initialization
Lou(k,j)=2; %initialization

while Lou(k,j)<L
 uu(k+1,j)=fsolve...
     .
     .

The problem is how to use fsolve in solving for the value of uu for each time it goes into the while loop; k is the counter

[jgg]

This isn't a minimal example. How is Lou(k,j) updated? What is the solving criterion? You're making this very difficult for anyone to help you with.

For instance, the code you've provided doesn't even run because uuf{k,j}(2) is a single value, not a vector. Additionally, you're not pre-declaring any of your variables, so they change size every iteration which makes it impossible to tell how anything should be structured.

[Isaac]

k is a counter for the while loop. I am solving for uu as long as the condition in the while loop is valid. So the length of k depends on the iteration steps.

Lou is updated in the while loop.

Lou(k+1,j)=Lou(k,j)+uu(k,j)

My

[Isaac]

I have tried using the fsolve this way also

testfunLu=@(x)eval(subs(funLu(k,j),uus,x))-p(k,j);

uuf(k,j)=fsolve(testfunLu,55);

[Isaac]

This is the pseudo code:

for j=1:N

      k=1;
      uu(k,j)=0;
      Lou(k,j)=2;
      p(k,j)=6000;
      pQ(k,j)=asin(1+Lou(k,j));
      Lh2(j)=200;
      k22(k,j)=Lh2(j)./Lou(k,j)

      if p(k,j)>pQ(k,j)
          while Lou(k,j)<Lh2
              k22f(k,j)=Lh2(j)./(Lou(k,j)+uus);
              funLu(k,j)=asin(1+k22f(k,j));
              testfunLu=@(x)eval(subs(funLu(k,j),uus,x))-p(k,j);
              uuf(k,j)=fsolve(testfunLu,55);
              uu(k+1,j)=uuf(k,j);
              Lou(k+1,j)=Lou(k,j)+uu(k+1,j);
              k22(k+1,j)=Lh2(j)./(Lou(k+1,j)+uu(k+1,j));
              p(k+1,j)=p(k,j);
              k=k+1
          end
      end
  end