Q = [1 0;
0 1/2];
R = 1/4
A = [0 1;
2 -1];
B = [0;
1];
tspan = [5 0];
P0 = 0;
P = ode45(@(t,P) odefun(t,P,Q,A,B,R), tspan, P0);
function dPdt = odefun(t,P,Q,A,B,R)
dPdt = -Q - A.'*P - P*A + P*B*R.^(-1)*B.'*P;
size(dPdt)
end
Q is 2 x 2, so addition or subtraction of Q and something else is going to give you at least a 2 x 2
A is 2 x 2, P is a scalar (because P0 is a scalar), so A.'*P is going to be 2 x 2.
P is a scalar, A is 2 x 2, so P*A is going to be 2 x 2
P is a scalar, B is 2 x 1 sp P*B is 2 x 1. R is a scalar, so P*B*R.^(-1) is 2 x 1. B is 2 x 1 so B.' is 1 x 2, and (2 x 1) * (1 x 2) is 2 x 2. P is scalar, so 2 x 2 * scalar is 2 x 2.
Thus you have multiple terms all of which are 2 x 2.
Your ode function must return a column vector, and the column vector must have the same number of elements as your P0 -- so MATLAB is expecting your function to return a scalar, not 2 x 2. Even if you reshape the 2 x 2 into a column vector, that would give you 4 elements, which would not match the single boundary conditon P0 that you have.
