Vectorizing equations for fsolve
Yes. You can simplify this code. I made two main changes. 

First, I forget that \$x_{50}\$ is known and write an equation for it just like all the other variables, then I replace that equation by the equation \$x_{50}=1\$.

Second, this is a system of linear equations. Notice that on both sides of the equation we have \$a_{ij}b_{kij}\$ with \$\{i,j,k\}=\{s,t,r\}\$ on the left and \$\{i,j,k\}=\{r,t,s\}\$ on the right, and both are summed over \$t\$ which corresponds to \$j\$. So first I find the matrix \$C_{ik}=\sum_ja_{ij}b_{kij}\$, which requires quite some tweaking due to the messed up dimensions/slices. Now the equation reads

\$\displaystyle x_{s}\sum_{r=1}^R C_{sr} = \sum_{r=1}^R C_{rs}x_r \$

for each \$s\in1,2,\ldots,S\$. The sum of the left is just sum(C,2) in Matlab, and the sum on the right is the matrix multiplication of \$C^t\$ with the vector of unknowns \$\vec x\$. We can write the equation in matrix form as

\$\displaystyle \vec c\vec x-C\vec x=0 \$

where \$[\vec c]_s=\sum_r C_{sr}\$. Written another way, we have

\$\displaystyle \left(\mathrm{diag}(\vec c)-C\right)\vec x=0. \$

Now since we know \$x_{50}=1\$ just replace the last row of the matrix with \$[0,0,\ldots,0,1]\$. Then if \$M=\mathrm{diag}(\vec c)-C\$ our equation is \$M\vec x=\vec v\$ where \$\vec v=[0,0,\ldots,0,1]^t\$. We can solve this system directly using \.

%% Define parameters
T = 3; % Number of types (T)
S = 3; % Number of senders (S)
R = S; % Number of receivers (R), equals the number of senders
A = [.1 .1 .1   % Columns are senders (S), rows are types (T)
     .3 .3 .6
     .6 .6 .3]; 
B1 = [.9 .0 .1     
      .5 .1 .2     
      .8 .2 .1];   
B2 = [.0 .7 .2     
      .2 .9 .0     
      .1 .8 .1];   
B3 = [.1 .3 .7     
      .3 .0 .8     
      .1 .0 .8];   

C = squeeze(sum(A.*cat(3,B1,B2,B3),1)).';

M = diag(sum(C,1))-C; % matrix equations
M(end,:) = [zeros(1,S-1) 1]; % replace the last equation

v = [zeros(S-1,1);1]; % modify final RHS entry

F = M\v;


An additional benefit is that solving linear systems is much faster than solving nonlinear systems, so this code is much faster. I couldn't test with different R, S, or T so something might have to change, I'm not sure if it will just work, but I think it should.



With the new method for generating the equations, we can make a few more simplifications to avoid having to concatenate the Bi's explicitly, using a cell array:

T = 3; % Number of types (T)
S = 3; % Number of senders (S)
R = S; % Number of receivers (R), equals the number of senders
A = [.1 .1 .1   % Columns are senders (S), rows are types (T)
     .3 .3 .6
     .6 .6 .3]; 
B{1} = [.9 .0 .1     
      .5 .1 .2     
      .8 .2 .1];   
B{2} = [.0 .7 .2     
      .2 .9 .0     
      .1 .8 .1];   
B{3} = [.1 .3 .7     
      .3 .0 .8     
      .1 .0 .8];   

C = squeeze(sum(A.*cat(3,B{:}),1)).';

M = diag(sum(C,1))-C;
M(end,:) = [zeros(1,S-1) 1];

v = [zeros(S-1,1);1];

F = M\v;