If you have symbolic toolbox, you can use solve()
syms M S R
Rv = rand(10, 1);
eq1 = M == S*3;
eq2 = S == R+M;
sol = solve([eq1; eq2], [M S]);
M_sol = subs(sol.M, Rv);
S_sol = subs(sol.S, Rv);
Result
>> M_sol
M_sol =
-14923794479723871/18014398509481984
-4255113867394491/4503599627370496
-108056045000625/2251799813685248
-16610538751653585/18014398509481984
-9792936758686299/18014398509481984
-1338449427099159/18014398509481984
-1653620415623073/2251799813685248
-5201938487820639/18014398509481984
-1662959758749849/9007199254740992
-694097603644977/2251799813685248
>> S_sol
S_sol =
-4974598159907957/18014398509481984
-1418371289131497/4503599627370496
-36018681666875/2251799813685248
-5536846250551195/18014398509481984
-3264312252895433/18014398509481984
-446149809033053/18014398509481984
-551206805207691/2251799813685248
-1733979495940213/18014398509481984
-554319919583283/9007199254740992
-231365867881659/2251799813685248
If you want floating-point answer
M_sol = double(M_sol);
S_sol = double(S_sol);
Reault
>> M_sol
M_sol =
-0.2198
-0.2836
-0.0640
-0.9528
-0.4228
-0.8079
-1.0427
-0.7487
-0.8037
-0.6678
>> S_sol
S_sol =
-0.0733
-0.0945
-0.0213
-0.3176
-0.1409
-0.2693
-0.3476
-0.2496
-0.2679
-0.2226
Alternatively use fsolve():
R = rand(10, 1);
eq = @(M,S) [M-S*3; S-(R+M)];
sol = fsolve(@(x) eq(x(1:10),x(11:20)), zeros(20, 1));
M_sol = sol(1:10);
S_sol = sol(11:20);
