solving two maximization non linear problems

Hello, I'm trying to solve to maximization problems: Utility_A and Utility_B, using the following code. Utility_A maximize by changing quantities: Epsilon_A_S,Epsilon_A_c. Utility_B maximize by changing quantities: Epsilon_B_S,Epsilon_B_c.

function constraints:

0<=Epsilon_A_S,Epsilon_B_S<=2
-1<=Epsilon_B_c,Epsilon_B_c<=1

Once these quantities sets by each utility function, an equilibrium rule sets:

if Epsilon_A_c=-Epsilon_B_c and Epsilon_A_S+Epsilon_B_S=2 stop the algorithm
% else...
if fitted_Epsilon_A_S+fitted_Epsilon_B_S>2
    S_0= S_0+0.00001;
else
    if fitted_Epsilon_A_S+fitted_Epsilon_B_S<2
        S_0= S_0-0.00001;
    end
end

Total_Epsilon_c=fitted_Epsilon_A_c+fitted_Epsilon_B_c;
if Total_Epsilon_c>0
    Call_Price_0= Call_Price_0+0.00001;
else
    if Total_Epsilon_c<0
        Call_Price_0= Call_Price_0-0.00001;
    end
end

and maximize two utility functions again until it converge. The output result should set Call_Price_0,S_0, Epsilon_B_S,Epsilon_B_c and Epsilon_A_S,Epsilon_A_c. The numerical should be: Epsilon_A_S=~2,Epsilon_B_S=~0, Epsilon_A_c=-1,Epsilon_B_c=1, S_0=~10.0 , Call_Price_0=~0.7. The Code does not converge and the parameters are not optimal.

Code Here:

clear all; close all
CP=[];
Cond3=0;
S=10;
H_mu=0.8;
S_A=S+H_mu;
S_B=S-H_mu;
Sigma_i=1.5;
K=10;
T=1;
r=0.02;
Gamma=0.1;
V=2;I=2;

S_0=(S-Gamma*(Sigma_i^2)*(V/I))/(1+r)
z=(K-S)/Sigma_i
Call_Price_0=((S-Gamma*(Sigma_i^2)*(V/I)-K)*(1-normcdf(z+Gamma*Sigma_i*(V/I)))+normpdf(z+Gamma*Sigma_i*(V/I))*Sigma_i)/(1+r)

while  Cond3<2 

S_Initial=S_0;Call_Price=Call_Price_0;
A=@(Epsilon_A_S,Epsilon_A_c) ...
    (((V/I)-Epsilon_A_S)*S_Initial-Epsilon_A_c*Call_Price)*(1+r)+Epsilon_A_S*S_A-(Gamma/2)*(Epsilon_A_S^2)*(Sigma_i^2);
B=@(Epsilon_A_S,Epsilon_A_c) ...
    (((V/I)-Epsilon_A_S)*S_Initial-Epsilon_A_c*Call_Price)*(1+r)+(Epsilon_A_S+Epsilon_A_c)*S_A-(Gamma/2)*((Epsilon_A_S+Epsilon_A_c)^2)*(Sigma_i^2)-Epsilon_A_c*K;
z_A=(K-S_A)/Sigma_i;
Utility_A=@(Epsilon_A_S,Epsilon_A_c)...
    -exp(-Gamma*A(Epsilon_A_S,Epsilon_A_c))*normcdf(z_A+Gamma*Epsilon_A_S*Sigma_i)+...
    -exp(-Gamma*B(Epsilon_A_S,Epsilon_A_c))*(1-normcdf(z_A+Gamma*(Epsilon_A_S+Epsilon_A_c)*Sigma_i));

%Take the negative of the utility. If we minimize that, we maximize
%utitliy. This function is unconstrained, but when either epsilon goes over
%the bound value, the cost function value goes to inf.
minimize_function=@(Epsilon_A_S,Epsilon_A_c)...
    -Utility_A(Epsilon_A_S,Epsilon_A_c)+...
    1/(Epsilon_A_S<=2)-1+...
     1/(Epsilon_A_S>=0)-1+...
    1/(Epsilon_A_c>=-1)-1;
intial_Epsilon_A_S=(V/I);intial_Epsilon_A_c=0;
fitted_val=...
    fminsearch(@(val)minimize_function(val(1),val(2)),...
    [intial_Epsilon_A_S,intial_Epsilon_A_c]);
[fitted_Epsilon_A_S,fitted_Epsilon_A_c]=deal(fitted_val(1),fitted_val(2));
%%Player B
C=@(Epsilon_B_S,Epsilon_B_c) ...
    (((V/I)-Epsilon_B_S)*S_Initial-Epsilon_B_c*Call_Price)*(1+r)+Epsilon_B_S*S_B-(Gamma/2)*(Epsilon_B_S^2)*(Sigma_i^2);
D=@(Epsilon_B_S,Epsilon_B_c) ...
    (((V/I)-Epsilon_B_S)*S_Initial-Epsilon_B_c*Call_Price)*(1+r)+(Epsilon_B_S+Epsilon_B_c)*S_B-(Gamma/2)*((Epsilon_B_S+Epsilon_B_c)^2)*(Sigma_i^2)-Epsilon_B_c*K;
z_B=(K-S_B)/Sigma_i;
Utility_B=@(Epsilon_B_S,Epsilon_B_c)...
    -exp(-Gamma*C(Epsilon_B_S,Epsilon_B_c))*normcdf(z_B+Gamma*Epsilon_B_S*Sigma_i)+...
    -exp(-Gamma*D(Epsilon_B_S,Epsilon_B_c))*(1-normcdf(z_B+Gamma*(Epsilon_B_S+Epsilon_B_c)*Sigma_i));

%Take the negative of the utility. If we minimize that, we maximize
%utitliy. This function is unconstrained, but when either epsilon goes over
%the bound value, the cost function value goes to inf.
minimize_functionB=@(Epsilon_B_S,Epsilon_B_c)...
    -Utility_B(Epsilon_B_S,Epsilon_B_c)+...
    1/(Epsilon_B_S<=2)-1+...
     1/(Epsilon_B_S>=0)-1+...
    1/(Epsilon_B_c>=-1)-1;
intial_Epsilon_B_S=(V/I);intial_Epsilon_B_c=0;
fitted_val=...
    fminsearch(@(val)minimize_functionB(val(1),val(2)),...
    [intial_Epsilon_B_S,intial_Epsilon_B_c]);
[fitted_Epsilon_B_S,fitted_Epsilon_B_c]=deal(fitted_val(1),fitted_val(2));

%%Equilibrium Conditions
if fitted_Epsilon_A_S+fitted_Epsilon_B_S>2
    Cond_Epsilpn_S=0;
    S_0= S_0+0.00001;
else
    if fitted_Epsilon_A_S+fitted_Epsilon_B_S<2
        Cond_Epsilpn_S=0;
        S_0= S_0-0.00001;
    else
            Cond_Epsilpn_S=1;
    end
end

Total_Epsilon_c=fitted_Epsilon_A_c+fitted_Epsilon_B_c;
if Total_Epsilon_c>0
    Cond_Epsilpn_c=0;
    Call_Price_0= Call_Price_0+0.00001;
else
    if Total_Epsilon_c<0
        Cond_Epsilpn_c=0;
        Call_Price_0= Call_Price_0-0.00001;
    else
            Cond_Epsilpn_c=1;
    end
end
Cond3=Cond_Epsilpn_c+Cond_Epsilpn_S
CP=[CP; S_0 Call_Price fitted_Epsilon_A_c]
end

[Yossi]

Please note that the convergence in this code does not meet any criteria for Equilibrium Conditions