MATLAB: Does solve() after sym(‘x’, ‘positive’) return a negative number

solve sym positive

I'm trying to solve a system of equations. I have defined 4 variables to be positive, but solve() returns a negative value. As far as I understand sym('x', 'positive') command, the result of solve should be an empty set if the only solution is negative… Am I missing something? Here is my code:
syms l_1 s_1 s_2 l_2 clear;
l_1 = sym('l_1', 'positive');
l_2 = sym('l_2', 'positive');
s_1 = sym('s_1', 'positive');
s_2 = sym('s_2', 'positive');
ft_1 = 10;
ft_2 = 10;
g_1 = 10;
beta = 0.9;
lambda = 0.6;
gama_0 = 0;
gama_1 = 0.5;
gama_2 = 0.2;
gama_3 = 0.2;
gama_4 = 0.2;
a = 1;
alfa_0 = 1;
alfa_1 = 0.5;
alfa_2 = 0.5;
L = 8;
delta_1 = 2;
delta_2 = 2;
V1 = log(ft_1 + g_1 + l_1) - delta_1*(l_1)^2 - delta_2*(s_1)^2+ beta*(log(ft_2 + lambda*(gama_0 + gama_1*s_1 + gama_2*l_1 + gama_3*a) + l_2) - delta_1*(l_2)^2 - delta_2*(s_2)^2) + beta^2*(log(L*alfa_0*((0.5*(gama_0 + gama_1*s_1 + gama_2*l_1 + gama_3*a) + 0.5*(gama_0 + gama_1*s_2 + gama_2*(l_1 + l_2) + gama_3*a + gama_4*(gama_0 + gama_1*s_1 + gama_2*l_1 + gama_3*a)))^alfa_1)*(l_1 + l_2)^alfa_2)- delta_1*(L)^2);
pVps_1 = diff(V1,s_1)
pVpl_1 = diff(V1,l_1)
pVps_2 = diff(V1,s_2)
pVpl_2 = diff(V1,l_2)
[xl1, xl2, xs1, xs2] = solve(pVpl_1, pVpl_2,pVps_1, pVps_2, l_1, l_2, s_1, s_2)
Results that I get are:
xl1 = -0.15635731765851100895734570889025
xl2 = -0.22318136524599875462767864215809
xs1 = 0.13647789503117958882560294058743
xs2 = 0.12006738235828457981821794454359

Best Answer

I recommend against using ls and 1s (ELs and ONEs) as your variable names since it's very hard to read.
The solver is returning negative answers because it cannot find an analytical solution. Once it cannot find an analytical solution it uses a numeric solver that does not obey the constraints.
You can explicitly call the numeric solver using:
evalin(symengine, 'numeric::fsolve([pVpl_1, pVpl_2,pVps_1, pVps_2],[l_1=0..infinity,l_2=0..infinity,s_1=0..infinity,s_2=0..infinity])')
But this is failing for me.