Troubles with interaction problem

Hey guys! I'm having trouble finding a value for "beta" in the following code:

 e=1.239e-7;
 c=5.154e-4;
 d=4.823;

 me=2.0475e8;         
 mv=28.96;              

 se=0.06.*me;            
 sv=6.081;              

 imax=15;                %Max number of interactions
 tol=0.001;              %Tolerance for convergence

    Edes(1)=me;
    Vdes(1)=mv;

    Ered_des(1)=(Edes(1)-me)./se;
    Vred_des(1)=(Vdes(1)-me)./se;

    dg_de(1)=se./(me+se.*Ered_des(1));
    dg_dv(1)=-2.*c.*sv.*(mv+sv.*Vred_des(1))./(d+c.*(mv+sv.*Vred_des(1)).^2);

    alfa_e(1)=(dg_de(1))./(sqrt(dg_de(1).^2+dg_dv(1).^2));
    alfa_v(1)=(dg_dv(1))./(sqrt(dg_de(1).^2+dg_dv(1).^2));

    syms beta
    g(1)=log( e.*((-alfa_e(1).*beta(1)).*se+me)./(c.*(-alfa_v(1).*beta(1).*sv+mv).^2+d) );
    [beta(1)] = solve(g(1)==0);
    beta(1)=vpa((subs(beta(1))));

    for ii=1:imax

    beta(ii+1)=0;

    if abs(beta(ii+1)-beta(ii))>0.001
    Ered_des(ii+1)=-alfa_e(ii).*beta(ii);
    Vred_des(ii+1)=-alfa_v(ii).*beta(ii);
    dg_de(ii+1)=se./(me+se.*Ered_des(ii+1));
    dg_dv(ii+1)=-2.*c.*sv.*(mv+sv.*Vred_des(ii+1))./(d+c.*(mv+sv.*Vred_des(ii+1)).^2);

    alfa_e(ii+1)=(dg_de(ii+1))./(sqrt(dg_de(ii+1).^2+dg_dv(ii+1).^2));
    alfa_v(ii+1)=(dg_dv(ii+1))./(sqrt(dg_de(ii+1).^2+dg_dv(ii+1).^2));

    syms beta
    g(ii+1)=log( e.*((-alfa_e(ii+1).*beta(ii+1)).*se+me)./(c.*(-alfa_v(ii+1).*beta(ii+1).*sv+mv).^2+d) );
    [beta(ii+1)] = solve(g(ii+1)==0);
    beta(ii+1)=vpa((subs(beta(ii+1))));

    else
        Beta=double(beta(ii+1));

    end

    if ii==imax
         disp('Exceeded max. number of interactions!')
    end

 end

When I try this, the following error appears:

"Error using mupadmex Error in MuPAD command: Index exceeds matrix dimensions.

Error in sym/subsref (line 1577) B = mupadmex('symobj::subsref',A.s,inds{:});

Error in InvarianteHL (line 62) g(ii+1)=log( e.*((-alfa_e(ii+1).*beta(ii+1)).*se+me)./(c.*(-alfa_v(ii+1).*beta(ii+1).*sv+mv).^2+d) );"

Please please help?

[bym]

this line

for ii=1:imax beta(ii+1)=0;

is improper syntax maybe there is a missing return and you intend

for ii=1:imax 
beta(ii+1)=0;

[raquel]

It is like the way you wrote last on my code, the way I pasted it here messed it up :)

Answer 1

[bym]

here is a partial answer, without using symbolic tool box...see my comment above and you should be able to adapt the code to what you want

clear
e=1.239e-7; c=5.154e-4; d=4.823;

me=2.0475e8; mv=28.96;

se=0.06.*me; sv=6.081;

imax=15; %Max number of interactions

tol=0.001; %Tolerance for convergence

      Edes(1)=me;
      Vdes(1)=mv;
      Ered_des(1)=(Edes(1)-me)./se;
      Vred_des(1)=(Vdes(1)-me)./se;
      dg_de(1)=se./(me+se.*Ered_des(1));
      dg_dv(1)=-2.*c.*sv.*(mv+sv.*Vred_des(1))./(d+c.*(mv+sv.*Vred_des(1)).^2);
      alfa_e(1)=(dg_de(1))./(sqrt(dg_de(1).^2+dg_dv(1).^2));
      alfa_v(1)=(dg_dv(1))./(sqrt(dg_de(1).^2+dg_dv(1).^2));
      %syms beta
      g=@(beta)log( e.*((-alfa_e(1).*beta).*se+me)./(c.*(-alfa_v(1).*beta.*sv+mv).^2+d) );
      %[beta(1)] = solve(g(1)==0);
      %beta(1)=vpa((subs(beta(1))));
      b(1) = fzero(g,10);

[raquel]

From the line:

   syms beta
   g(ii+1)..

I replaced with:

     g=@(beta)log( e.*((-alfa_e(ii+1).*beta).*se+me)./(c.*(-alfa_v(ii+1).*beta.*sv+mv).^2+d) );
    beta(ii+1) = fzero(g,10);

    dif=abs(beta(ii+1)-beta(ii));

         if dif<0.001
        Beta=double(beta(ii+1))
         break
        end

And it works! :D Thank you.

But in this case I'd have to make a guess for beta, right? As you did using fzero(g,10).. Isn't there any other way to solve for beta?