MATLAB how do I access a specific coefficient in a symbolic equation system solution?

Question

I need to run a simple Monte Carlo varying coefficients on a system of equations. I need to record the solved coefficient of one of the variables each time.

The following gets me results from a single run:

syms alpha gamma Ps Pc beta lambda Pp Sp Ss Dp Ds;

eq1 = -Ss + alpha + 0.17*Ps - 1*Pc;
eq2 = -Sp + beta + 0.2*Pp;
eq3 = -Ds + gamma - 0.2*Ps + 1*Pp;
eq4 = -Dp + lambda - 0.17*Pp + 1*Ps;
eq5 = Ss - Ds;
eq6 = Sp - Dp;

ans1 = solve(eq1,eq2,eq3,eq4,eq5,eq6,'Ps','Pp','Ss','Ds','Sp','Dp');

disp('Ps')
vpa(ans1.Ps,3)
disp('Pp')
vpa(ans1.Pp,3)
disp('Ss')
vpa(ans1.Ss,3)
disp('Ds')
vpa(ans1.Ds,3)
disp('Sp')
vpa(ans1.Sp,3)
disp('Dp')
vpa(ans1.Dp,3)

I will be varying several of the variables (on Ps, Pp, and Pc), and recording the coefficient on Pc in each of the reduced form equations (i.e. the coefficient on Pc that shows up after vpa(ans1.xx)--so in the case above, it would be a 1x6 vector [-0.429,-1.16,-1.07,-1.07,-0.232,-0.429,-1.16]).

I'm very new to MATLAB, but I'm sure I can figure out how to implement the loop code to do the model iterations. What I can't figure out is how to record the vector of coefficients after each iteration. Is there some "accessor" that will give me just the one coefficient for each equation each time? Something like vpa(ans1.ps.coef(pc)) (which is a total shot in the dark, and it's wrong, but hopefully you get the idea).

Answer 1

There's probably a better way to do this, but this is all I could come up to at this moment.

Step 1: In order to obtain the coefficient of Pc as a double from ans1.Ps, you can use the subs function, as follows:

 subs(ans1.Ps,{alpha,Pc,beta,gamma,lambda},{0,1,0,0,0});

Step 2a:

To get a vector of all coefficients per ans1 expression (say ans1.Ps) you can use something like this:

N=numel(symvar(ans1.Ps)); % obtain number of coefs
cp=num2cell(eye(N));      % create a cell array using unit matrix, so each iteration a different coef will be selected

for n=1:N;
   coefs(n)=subs(ans1.Ps,{alpha,Pc,beta,gamma,lambda},cp(n,:));
end

Step 2b:

Alternatively, you want to get just Pc, but from all the ans1 expressions. If so then you can do the following:

SNames = fieldnames(ans1); % get names of ans1 expressions
for n = 1:numel(SNames) 
    expr = ans1.(SNames{n}); % get the expression itself
    pc(n)=subs(expr,{alpha,Pc,beta,gamma,lambda},{0,1,0,0,0}); % obtain just pc
end

You can now combine the two if you want all info about the coefficients.

Edit:

To store the retrieved Pc per iteration you can do the following:

alpha=[3 1 4 6 7] % just a vector of values
beta = [6 7 8 5 2]
SNames = fieldnames(ans1); % get names of ans1 expressions

for n = 1:numel(SNames) 
    expr = ans1.(SNames{n}); % get the expression itself

    for n1=1:numel(alpha)
        for n2=1:numel(beta) 

            pc(n,n1,n2)=subs(expr,{alpha,Pc,beta,gamma,lambda},{alpha(n1),1,beta(n2),0,0})

         end
    end
end