function computeEquilibrium_threePlayers
% Computes equilibrium strategie in a first-price auction with three
% players using the boundary value method with fixed-point iterations
%
%! When using this code please cite the following article:
%! G. Fibich and N. Gavish, Numerical simulations of asymmetric first-price auction, GEB, 2011
%! doi:10.1016/j.geb.2011.02.01

%% Define distribution functions
F1overf1=  @(x)(x) ;                                                        % Corresponding to power law distribution, F1=x
F2overf2=  @(x)(2*(exp(x/2)-1)) ;                                    % Corresponding to Weibull distribution, F2=c (1-exp(-x/2))
F3overf3=  @(x)(exp(x.^2).*erf(x)*sqrt(pi)/2);                 % Corresponding to truncated normal distribution, F3=c erf(x) 

%% Compute solution to (1) using  the fixed point iterations 
fixedPointIterations(F1overf1,F2overf2,F3overf3);
function fixedPointIterations(F1overf1,F2overf2,F3overf3)
%% Construct grid for v2
h=1/2000;
v3=(h:h:1-h)';
nGrid=numel(v3);

%% Construct finite-difference approximation for the differentiation operators
e=ones(nGrid,1);D = spdiags([-e 0*e e], -1:1, nGrid, nGrid)/2/h;
% Handle boundary v(1)=1 for Dv
Dv=D;Dv_RHS=e*0;
Dv_RHS(end)=1/2/h;
% Handle one-sided boundary b(1)=? for Db
Db=D;Db_RHS=e*0;
Db(end,end-2:end)=[1 -4 3]/2/h;

%% Construct initial guess
b=0.9*v3;
v1=v3;
v2=v3;

%% Compute solution using iterations
for iter=1:25
    % Calculate v1^{(k+1)}(v3)
    v_coeff=-F1overf1(v1)./F3overf3(v3)./((v3-b).*(v1+v2-2*b)-(v1-b).*(v2-b));
    LHS_v=Dv+spdiags(v_coeff.*(v2+v3-2*b),0,nGrid,nGrid);
    RHS_v=((v3+v2-2*b).*b+(v2-b).*(v3-b)).*v_coeff-Dv_RHS;
    v1=LHS_v\RHS_v;
    
    % Calculate v2^{(k+1)}(v3)
    v_coeff=-F2overf2(v2)./F3overf3(v3)./((v3-b).*(v1+v2-2*b)-(v1-b).*(v2-b));
    LHS_v=Dv+spdiags(v_coeff.*(v1+v3-2*b),0,nGrid,nGrid);
    RHS_v=((v3+v1-2*b).*b+(v1-b).*(v3-b)).*v_coeff-Dv_RHS;
    v2=LHS_v\RHS_v;

    % Calculate b^{(k+1)}(v3)
    b_coeff=2./F3overf3(v3).*(v1-b).*(v2-b)./((v1+v2-2*b).*(v3-b)-(v1-b).*(v2-b));
    LHS_b=Db+spdiags(b_coeff,0,nGrid,nGrid);
    RHS_b=v3.*b_coeff-Db_RHS;
    b=LHS_b\RHS_b;
    
end

%% Plot result
plot(b,v1,b,v2,b,v3);
box on;
legend('v_1','v_2','v_3','location','northwest');
xlabel('b');ylabel('v_i','rotation',0);
shg

%% Compute b_bar by extrapolation 
b_bar=[-1 4 -6 4]*b(end-3:end)
return