MATLAB code for the computations in "Nonlinear Pricing, market coverage, and competition" 
by Huanxing Yang and Lixin Ye (Theoretical Economics 3 (2008)), March 4, 2008


%This program produces the lowest type covered as k varies (uniform
%distribution)

%left end of time region
global k

k=1.3; %right end of time region
t0=0.5+sqrt(3*k)/4; %middle 
t0=t0+0.001
while k>=0.1
t0=t0-0.0001

error=1e-3; %suppose the error is within 1e-3, could change it.

[T,Y] = ode45(@vdp4,[t0, 1], [k/4, sqrt(3*k)/2]);

n=length(Y(:,2)); %the length of vector of y'
n

Y(n,2)

while (Y(n,2)>1.0+error | Y(n,2)<1.0-error)
    
    t0=t0-0.0001;

    [T,Y] = ode45(@vdp4,[t0, 1], [k/4, sqrt(3*k)/2]);

    n=length(Y(:,2)); %the length of vector of y'
    k
    t0 % initial time
    Y(n,2) % y' at time=1
    

end %while
t1=t0-sqrt(k/3);
plot(k,t1,'.');
hold on;
k=k-0.01
end
k=0.10:0.001:1.3;
y4=0.5-sqrt(3*k)/12;
plot (k,y4,'.');


%vdp4.m:

%differention equations for the competition range (uniform distribution)

function dy = rigid(t,y)
global k

dy = zeros(2,1);    % a column vector

dy(1)=y(2);
dy(2)=(y(2)^2-2*t*y(2)+2*y(1))/k+2;



%This program produces the lowest type covered as k varies (exponential
%distribution)

global k

k=1.3;
t0=0.355;
t1=0.975;
error=1e-2;
while k>=0.1


[T,Y] = ode45(@vdp5, [t0, t1], [0.000001, 0]);

n=length(Y(:,2));
Y(n,2)
t1=t1-0.05
while (Y(n,2)>1.0+t1-exp(1-t1)+error | Y(n,2)<1.0+t1-exp(1-t1)-error)
    t0=t0+0.001;
    t1=t1-0.05
    [T,Y] = ode45(@vdp5, [t0, t1], [0.000001, 0]);
    n=length(Y(:,2))
    Y(n,1)
    while (Y(n,1)<k/4-error | Y(n,1)>k/4+error)
    t1=t1+0.001;
    
    [T,Y] = ode45(@vdp5, [t0, t1], [0.000001, 0]);
    n=length(Y(:,2));
    
    k
    t0
    t1
    Y(n,2)
    Y(n,1)
end
end
plot(k,t0, '.');
hold on;
k=k-0.01;
end

hold on;

k=1.3;
t0=0.365;
t1=0.95
error=1e-2;
while k>=0.1
if k>=0.4
    t0=t0-0.01;
else
    t0=t0-0.015;
end
t1=t1-0.03;


[T,Y] = ode45(@vdp5, [t0, t1], [0.000001, 0]);

n=length(Y(:,1));
m1=Y(n,1);
m2=Y(n,2);

[T,Y] = ode45(@vdp6, [t1, 1], [m1, m2]);
n=length(Y(:,2));

while (Y(n,2)>1+error | Y(n,2)<1-error)
    t0=t0+0.001;
    t1=t1-0.03;
    
    [T,Y] = ode45(@vdp5, [t0, t1], [0.000001, 0]);
    n=length(Y(:,2));
    
    while (Y(n,1)<k/4-error | Y(n,1)>k/4+error)
    t1=t1+0.001;
    
    [T,Y] = ode45(@vdp5, [t0, t1], [0.000001, 0]);
    n=length(Y(:,2));
    k
    t0
    t1
    m1=Y(n,1)
    m2=Y(n,2)
end

[T,Y] = ode45(@vdp6, [t1, 1], [m1, m2]);
n=length(Y(:,2));
end
plot(k,t0, '.');
hold on;
k=k-0.01;
end



%this program produces quality schedules (exponential distribution k=0.4)

global k

k=0.4;
t0=0.451;
t1=t0+0.35;
error=1e-2;

[T,Y] = ode45(@vdp5, [t0, t1], [0.000001, 0]);

n=length(Y(:,1));
Y(n,1)
    while (Y(n,1)<k/4-error | Y(n,1)>k/4+error)
    t1=t1+0.0001;
    
    [T,Y] = ode45(@vdp5, [t0, t1], [0.000001, 0]);
    n=length(Y(:,2));
    
    t0
    t1
    Y(n,1)
    Y(n,2)
end
z=1+t1-exp(1-t1)-Y(n,2);
z
t0
t1
plot(T,Y(:,2));
hold on;
x=t1:0.005:1;
y=1+x-exp(1-x);
plot(x,y); 
hold on;

t0=0.326;
t1=t0+0.3;


[T,Y] = ode45(@vdp5, [t0, t1], [0.000001, 0]);

n=length(Y(:,1));
Y(n,1)
    while (Y(n,1)<k/4-error | Y(n,1)>k/4+error)
    t1=t1+0.0001;
    
    [T,Y] = ode45(@vdp5, [t0, t1], [0.000001, 0]);
    n=length(Y(:,2));
    
    t0
    t1
    m1=Y(n,1)
    m2=Y(n,2)
end
t1
plot(T,Y(:,2));
hold on;

[T,Y] = ode45(@vdp6, [t1, 1], [m1, m2]);
n=length(Y(:,2));
z=1-Y(n,2);
z
plot(T,Y(:,2)); 



%vdp5.m

%differention equations for the local monopoly range (exponential distribution)

function dy = rigid(t,y)
global k

dy = zeros(2,1);

dy(1)=y(2);
dy(2)=3-(1/2)*y(2)^2/(y(1))+t-y(2);


%vdp6.m

%differention equations for the competition range (exponential
%distribution)

function dy = rigid(t,y)
global k

dy = zeros(2,1);    % a column vector

dy(1)=y(2);
dy(2)=(y(2)^2-2*t*y(2)+2*y(1))/k+2+t-y(2);




%This program produces the lowest type covered as k varies (generalized
%uniform distribution)

global k

k=1.3;
t0=0.355;
t1=0.975;
error=1e-2;
while k>=0.1


[T,Y] = ode45(@vdp7, [t0, t1], [0.000001, 0]);

n=length(Y(:,2));
Y(n,2)
t1=t1-0.05
while (Y(n,2)>t1-(1-t1^2)/(2*t1)+error | Y(n,2)<t1-(1-t1^2)/(2*t1)-error)
    t0=t0+0.001;
    t1=t1-0.05
    [T,Y] = ode45(@vdp7, [t0, t1], [0.000001, 0]);
    n=length(Y(:,2))
    Y(n,1)
    while (Y(n,1)<k/4-error | Y(n,1)>k/4+error)
    t1=t1+0.001;
    
    [T,Y] = ode45(@vdp7, [t0, t1], [0.000001, 0]);
    n=length(Y(:,2));
    
    k
    t0
    t1
    Y(n,2)
    Y(n,1)
end
end
plot(k,t0, '.');
hold on;
k=k-0.01;
end

hold on;

k=1.3;
t0=0.365;
t1=0.96
error=1e-2;
while k>=0.1
if k>=0.4
    t0=t0-0.01;
else
    t0=t0-0.015;
end
t1=t1-0.03;


[T,Y] = ode45(@vdp7, [t0, t1], [0.000001, 0]);

n=length(Y(:,1));
m1=Y(n,1);
m2=Y(n,2);

[T,Y] = ode45(@vdp8, [t1, 1], [m1, m2]);
n=length(Y(:,2));

while (Y(n,2)>1+error | Y(n,2)<1-error)
    t0=t0+0.001;
    t1=t1-0.03;
    
    [T,Y] = ode45(@vdp7, [t0, t1], [0.000001, 0]);
    n=length(Y(:,2));
    
    while (Y(n,1)<k/4-error | Y(n,1)>k/4+error)
    t1=t1+0.001;
    
    [T,Y] = ode45(@vdp7, [t0, t1], [0.000001, 0]);
    n=length(Y(:,2));
    k
    t0
    t1
    m1=Y(n,1)
    m2=Y(n,2)
end

[T,Y] = ode45(@vdp8, [t1, 1], [m1, m2]);
n=length(Y(:,2));
end
plot(k,t0, '.');
hold on;
k=k-0.01;
end



%this program produces quality schedules (generalized uniform distribution k=0.4)


global k

k=0.3;
t0=0.482;
t1=t0+0.3;
error=1e-2;

[T,Y] = ode45(@vdp7, [t0, t1], [0.000001, 0]);

n=length(Y(:,1));
Y(n,1)
    while (Y(n,1)<k/4-error | Y(n,1)>k/4+error)
    t1=t1+0.0001;
    
    [T,Y] = ode45(@vdp7, [t0, t1], [0.000001, 0]);
    n=length(Y(:,2));
    
    t0
    t1
    Y(n,1)
    Y(n,2)
end

t0
t1
plot(T,Y(:,2));
hold on;
x=t1-0.005;
while x<=1-0.005
x=x+0.001
y=x-(1-x^2)/(2*x);
plot(x, y); 
hold on;
end


t0=0.314;
t1=t0+0.26;


[T,Y] = ode45(@vdp7, [t0, t1], [0.000001, 0]);

n=length(Y(:,1));
Y(n,1)
    while (Y(n,1)<k/4-error | Y(n,1)>k/4+error)
    t1=t1+0.0001;
    
    [T,Y] = ode45(@vdp7, [t0, t1], [0.000001, 0]);
    n=length(Y(:,2));
    t0
    t1
    m1=Y(n,1)
    m2=Y(n,2)    
end



%vdp7.m

%differention equations for the local monopoly range (generalized uniform
%distribution)

function dy = rigid(t,y)
global k

dy = zeros(2,1);

dy(1)=y(2);
dy(2)=3-(1/2)*y(2)^2/(y(1))+(t-y(2))/t;



%vdp8.m

%differention equations for the competition range (generalized uniform
%distribution)

function dy = rigid(t,y)
global k

dy = zeros(2,1);    % a column vector

dy(1)=y(2);
dy(2)=(y(2)^2-2*t*y(2)+2*y(1))/k+2+(t-y(2))/t;