Problem with matrix dimensions

Hello, my code is giving me the error :

Error using \

Matrix dimensions must agree.

I don't know why because the matrices I am trying to divide have the dimensions of 5x1 and 5x5

F = @sistema;
JF = @jacobiano;
y0 = zeros(5,1);
e = 1e-5;
nmax = 30;
[x, it, norma_diferenca] = metodo_newton2(F, JF, y0, e, nmax) 

x = 5x1

1.0e+04 * 0.3333 0 0 0 9.6667

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

it = 2

norma_diferenca = 0

function F = sistema(X)
beta1 = 0.0000009;
beta2 = 0.00000093;
M = 1;
p = 0.1;
mu = 0.00001;
rho = 0.00002;
theta = 0.005;
omega = 0.05;
gamma = 0.001;
alpha1 = 0.0001;
alpha2 = 0.004;
S = X(1);
E = X(2);
I = X(3);
Q = X(4);
V = X(5);
F(1) = p*M - beta1*S*E - beta2*S*I - (rho + mu)*S;
F(2) = beta1*S*E + beta2*S*I - (omega + alpha1 + mu)*E;
F(3) = omega*E - (gamma + alpha2 + mu)*I;
F(4) = alpha1*E + alpha2*I - (theta + mu)*Q;
F(5) = (1 - p)*M + rho*S + gamma*I + theta*Q - mu*V;
F = F(:);
end
function JF = jacobiano(X)
beta1 = 0.0000009;
beta2 = 0.00000093;
M = 1;
p = 0.1;
mu = 0.00001;
rho = 0.00002;
theta = 0.005;
omega = 0.05;
gamma = 0.001;
alpha1 = 0.0001;
alpha2 = 0.004;
S = X(1);
E = X(2);
I = X(3);
Q = X(4);
V = X(5);
JF = [-beta1*E - beta2*I - rho - mu, -beta1*S, -beta2*S, 0, 0;
    beta1*E, beta1*S - (omega + alpha1 + mu), beta2*S, 0, 0;
    0, omega, -(gamma + alpha2 + mu), 0, 0;
    0, alpha1, alpha2, -(theta + mu), 0;
    rho, 0, gamma, theta, -mu];
end

These are the codes of the two functions that qive the matrices

function [x, it, norma_diferenca] = metodo_newton2(F, JF, y0, e, nmax)
x = y0;
it = 0;
norma_diferenca = Inf;
while norma_diferenca > e && it < nmax
    Fx = F(x);
    JFx = JF(x);
    div = JFx \ (-Fx);
    x = x + div;
    norma_diferenca = norm(div);
    it = it + 1;
end
if it >= nmax
    disp('Atingido o número máximo de iterações sem convergência.');
end
end

And this is the code taht divide the two matrices.

Answer 1

[Torsten]

In "sistema", I made F a column vector by adding a last line as

F = F(:);

and it works.

Answer 2

[John D'Errico]

I'm not sure why you would use a Newton method to solve what solve will do directly, in an apparent algebraic form.

beta1 = 0.0000009;
beta2 = 0.00000093;
M = 1;
p = 0.1;
mu = 0.00001;
rho = 0.00002;
theta = 0.005;
omega = 0.05;
gamma = 0.001;
alpha1 = 0.0001;
alpha2 = 0.004;
syms S E I Q V
F(1) = p*M - beta1*S*E - beta2*S*I - (rho + mu)*S;
F(2) = beta1*S*E + beta2*S*I - (omega + alpha1 + mu)*E;
F(3) = omega*E - (gamma + alpha2 + mu)*I;
F(4) = alpha1*E + alpha2*I - (theta + mu)*Q;
F(5) = (1 - p)*M + rho*S + gamma*I + theta*Q - mu*V;

sol = solve(F,[S E I Q V])

sol = struct with fields:
    S: [2x1 sym]
    E: [2x1 sym]
    I: [2x1 sym]
    Q: [2x1 sym]
    V: [2x1 sym]

[sol.S,sol.E,sol.I,sol.Q,sol.V]

ans = 

The first solution has three of the unknowns as identically zero, with only S and V non-zero. That is probably a degenerate solution that you don't care to find. But the second solution is non-trivial.

If I look at the equations carefully, It seems there are only two equations that are not purely linear in the parameter. The first two. And there we see you have the terms S*E and S*I as products.

F(1) = p*M - beta1*S*E - beta2*S*I - (rho + mu)*S;
F(2) = beta1*S*E + beta2*S*I - (omega + alpha1 + mu)*E;

But, if we just add those first two equations, so form F(1) + F(2), the products go completely away. My guess is this makes the problem now almost solvable by hand. In the end, it will reduce to a quadratic, and so, you end up with two roots.