How can I solve a nonlinear matrix equation, and get the possible _matrices_ as my output, instead of a struct of 1-dimensional solutions?

I have searched online vigorously for answers to this question, and I couldn't find anything helpful. The closest I found is an unanswered question from 2019, at: Solve for an unknown matrix like solve() for dimension one - MATLAB Answers - MATLAB Central (mathworks.com)

I have a matrix equation, where all is known except for one matrix, P in this case. This equation has two possible solutions.

I can easily solve for the entries of P, and re-assemble them to different matrices P1 and P2, see below:

Here is the same code so you don't have to enter it manually:

clearvars; clc;

% setup
syms rho p1 p2 p4;
A = [0 1; 0 0];
B = [0;1];
H = [0 1];
R = rho;
Q = 1;
P = [p1 p2; p2 p4]; % P is symmetric

% solve this equation
eqns = A.'*P + P*A + H.'*Q*H - P*B * (R \ B.') * P == zeros(2,2);
[p1 p2 p4] = solve(eqns,[p1 p2 p4])

% Now I have to separate the results manually
P1(1,1) = p1(1);
P1(1,2) = p2(1);
P1(2,1) = p2(1);
P1(2,2) = p4(1)

P2(1,1) = p1(2);
P2(1,2) = p2(2);
P2(2,1) = p2(2);
P2(2,2) = p4(2)

However, are there no better ways to do this? I would like to solve for the matrix P directly instead. I would like the result to come out in matrix form. Since we have two possible solutions, the answer can be a structure or a three-dimensional matrix. That's not an issue as long as it's arranged in matrix form already.

Is this possible to do?

To be clear, I'm not looking for direct-form answers to my equation. I'm not looking for a (oh, flip that A around, multiply that B by an inverse, factorize the P and you get it in terms of all the other known matrices). I'm looking for a general form for any matrix equation, or, better yet, system of matrix equations.

[Torsten]

I would like to solve for the matrix P directly instead. I would like the result to come out in matrix form.

Already the Rolling Stones knew: You can't always get what you want. And this is also true here.

Answer 1

[Matt J]

I would like to solve for the matrix P directly instead ... Is this possible to do?

No, but the manual work of assembling P at the end can be reduced considerably, e.g., as below:

syms rho ; syms p [1,3];
A = [0 1; 0 0];
B = [0;1];
H = [0 1];
R = rho;
Q = 1;
P = [p1 p2; p2 p3]; % P is symmetric

% solve this equation
eqns = A.'*P + P*A + H.'*Q*H - P*B * (R \ B.') * P == zeros(2,2);
[pc{1:3}]=solve(eqns,p);


P=[pc{:,[1,2,2,3]}];

P1= reshape( P(1,:),2,2)

P1 = 

P2= reshape( P(2,:)  ,2,2)

P2 = 

[Ali]

I am blown away. I never thought to access elements of an array this way:

 X = [ 1 2 3]
X =
     1     2     3
>> X(:,[1 1 3 2 2])
ans =
     1     1     3     2     2

Thanks for the new knowledge! :-)