Problem Statement:
The problem is to find p and t1 functions with the right properties.
We want them both increasing in m1,m2. They cannot depend on
v1 and v2, except indirectly, through the players’ choice of m1,m2.
The probability p should be bounded between 0 and 1, and the tax t1
between 0 and 100. Also, we want a unique solution to exist when the
players maximize their payoffs and we want the equilibrium probability
to be arbitrarily close to 1. And our goal is to get t1 as close as possible
to v1/(v1+v2). We can limit v1 and v2 to be in [0,100] if necessary.

Update: Here we don’t want to SOLVE for p (we should make that clear) we just want to chose p qualitatively and then find m1^* and m2^*.

Eric: the PDE solutions below are exceptional cases where the m1^* is free and only m2^* is determined...these are special exceptional functions.

Note we cannot solve (in complete generality) for v1 and v2 in terms of t1 and t2 since t2 = 1 - t1, so only one parameter in general.

New version (we chose p[m1,m2] as a given function and simply solve for m1,m2).

Example 1

p[m1_,m2_]=1-b/(m1 m2)^b
t1[m1_,m2_]=m1/(m1+m2)

Here we are first treating v1 and v2 as constant in m1 and m2.

Eqs={D[p[m1,m2],m1](v1-t1[m1,m2])-p[m1,m2]D[t1[m1,m2],m1],D[p[m1,m2],m2](v2-1+t1[m1,m2])-p[m1,m2]D[1-t1[m1,m2],m2]}// FullSimplify

Solve[Eqs==0,{m1,m2}]

Eqs2={(b m1 m2-(((m1 m2))^(b+1))+(b^2) (m1+m2) (m1 (-1+v1)+m2 v1)), (b m1 m2- (((m1 m2))^(b+1))+(b^2) (m1+m2) (m2 (-1+v2)+m1 v2))}

Solve[Eqs2==0,{m1,m2}]

x=m1 m2  and y=m1+m2

xytom=Solve[{x==m1 m2,y==m1+m2},{m1,m2}]//FullSimplify

Note these are equivalent solutions (just switched order)

Eqs3={b x-x^(b+1)+b^2 y(-(1/2) (y+(Sqrt[-4 x+(y^2)]))+y v1),b x-x^(b+1)+b^2 y(-(1/2) (y-(Sqrt[-4 x+(y^2)]))+y v2)}

Eqs3p={b x-x^(b+1)+b^2 y(-(1/2) (y-(Sqrt[-4 x+(y^2)]))+y v1),b x-x^(b+1)+b^2 y(-(1/2) (y+(Sqrt[-4 x+(y^2)]))+y v2)}

Solve[Eqs3==0,{x,y}]

Reduce[Eqs3==0&& b>0,{x,y}]

Subtracting both equations and dividing by b^2 y we have...

repx=Solve[v1 y-(1/2) (y+(Sqrt[-4 x+(y^2)]))-(v2 y+(1/2) (-y+(Sqrt[-4 x+(y^2)])))==0,x][[1]]

Solve[v1 y-(1/2) (y+(Sqrt[-4 x+(y^2)]))-(v2 y+(1/2) (-y+(Sqrt[-4 x+(y^2)])))==0,y]
repy=%[[2]]

test={x-> .53,y->.331,v1->.76,v2->.42,b->1.3}
test2={x-> .53,y->.331,v1->1.76,v2->.42,b->1.3}

eqy0=b x-(x^(1+b))+(b^2) y (v1 y+(1/2) (-y-(Sqrt[-4 x+(y^2)])))/.repx//FullSimplify

eqx0=b x-(x^(1+b))+(b^2) y (v1 y+(1/2) (-y-(Sqrt[-4 x+(y^2)])))/.repy//FullSimplify

We reduce the above under the assumptions that v1>= v2 (Just so sqrt[(v1-v2)^2]=v1-v2 ow =v2-v1) and (v1-v2)<1 (just for single solutions in power function). Note that one square root needs to be <0 to avoid a 3

eqy1=(1/4) b (1-(((v1-v2))^2)) (y^2)-(4^(-1-b)) ((((1-(((v1-v2))^2)) (y^2)))^(1+b))-(1/2) (b^2) y (y-2 v1 y+(v1-v2) y)//FullSimplify
eqx1=b x-(x^(1+b))-((2 (b^2) x((1-2 v1)+(v1-v2)))/(1-(((v1-v2))^2)))//FullSimplify

Note the (1-3 v1+v2) should be (1- v1-v2) if v1>1+v2. Here is the check:

{eqx0,eqy0}/.#&/@{test,test2}
{eqx1,eqy1}/.#&/@{test,test2}

If (v1-v2)>1, then
eqx=
eqy=

eqy2=(1/4) b (1-(((v1-v2))^2))+(1/2) (b^2) (-1+v1+v2)-(4^(-1-b)) (((1-(((v1-v2))^2)))^(1+b))y^(2b)
eqx2=b-((2 (b^2) (1-v1-v2))/(1-(((v1-v2))^2)))-(x^b)

{eqx0,eqy0}/.#&/@{test,test2}
{eqx2 x,eqy2 y^2}/.#&/@{test,test2}

solx=Solve[eqx2==0,x][[1]]//FullSimplify

soly=Solve[eqy2==0,y][[1]]//FullSimplify

strsol={solx,soly};
strsol//MatrixForm

repsol={repx/.soly ,repy/.solx}//FullSimplify;
repsol//MatrixForm

solsxy={x->  (((b (1-(((v1-v2))^2)+2 b (-1+v1+v2))/(1-(((v1-v2))^2))))^(1/b)),y->2(((b (1-(((v1-v2))^2)+2 b (-1+v1+v2))/(1-(((v1-v2))^2))))^(1/(2b)))/(Sqrt[1-(((v1-v2))^2) ])}
solsxy1={x-> (((b-((2 (b^2)(1-v1-v2))/(1-(v1-v2)^2))))^(1/b)),y-> ((2 (((b-((2 (b^2)(1-v1-v2))/(1-(v1-v2)^2))))^(1/(2b))))/(Sqrt[1-(((v1-v2))^2) ]))}

Verification...

test3={v1->0.76`,v2->0.42`,b->1.3`};

#/.test3& /@ {{eqx0,eqy0},strsol,repsol,solsxy,solsxy1}

(Eqs3/.Flatten[#])/.test3 &/@{strsol,repsol,solsxy,solsxy1}//FullSimplify
(Eqs3p/.Flatten[#])/.test3 &/@{strsol,repsol,solsxy,solsxy1}//FullSimplify

xytom[[2]]/.solsxy//FullSimplify

msols={m1->(((1+(v1-v2))(((b-((2 (b^2)(1-v1-v2))/(1-(v1-v2)^2))))^(1/(2b))))/(Sqrt[1-(((v1-v2))^2)])),m2->(((1-(v1-v2))(((b-((2 (b^2)(1-v1-v2))/(1-(v1-v2)^2))))^(1/(2b))))/(Sqrt[1-(((v1-v2))^2)]))}

v20=.2;b0=8;
Plot[Evaluate[{(((((b-((2 (b^2) (1-v1-v2))/(1-(((v1-v2))^2)))))^((1/2)/b)) (1+v1-v2))/(Sqrt[1-(((v1-v2))^2)])),(((1-(v1-v2))(((b-((2 (b^2)(1-v1-v2))/(1-(v1-v2)^2))))^(1/(2b))))/(Sqrt[1-(((v1-v2))^2)]))}/.{b->b0,v2->v20}],{v1,Max[v20,1-v20]-.2,1+v20-.01},PlotRange->All]

Plot[(p[m1,m2]/.msols)/.{b->b0,v2->v20},{v1,Max[v20,1-v20]-.01,1+v20+.1},PlotRange->All]

Plot[(t1[m1,m2]/.msols)/.{b->b0,v2->v20},{v1,Max[v20,1-v20],1+v20-.01},PlotRange->All]

(p[m1,m2]/.msols)//FullSimplify//PowerExpand

(t1[m1,m2]/.msols)//FullSimplify

Solve[-1+(((v1-v2))^2)-2 b (-1+v1+v2)==0,v1]//FullSimplify

Reduce[0≤ 1-((-1+(((v1-v2))^2))/(-1+(((v1-v2))^2)-2 b (-1+v1+v2)))≤ 1&&b>0&&v1>v2>0]

Solving the first order equations for p. Here v1 and v2 and t1 can depend on m1,m2

DSolve[D[p[m1,m2],m1](v1[m1,m2]-t1[m1,m2])-p[m1,m2]D[t1[m1,m2],m1]==0,{p[m1,m2]},m1]

DSolve[D[p[m1,m2],m2](v2[m1,m2]-1+t1[m1,m2])-p[m1,m2]D[1-t1[m1,m2],m2]==0,{p[m1,m2]},m2]

compcond=(D[Log[(E^(InactiveIntegrate[-((Derivative[1,0][t1][K[1],m2])/(t1[K[1],m2]-v1[K[1],m2])), K[1], 1, m1])) (C[1])[m2]],m1,m2]-D[Log[(E^(InactiveIntegrate[-((Derivative[0,1][t1][m1,K[1]])/(-1+t1[m1,K[1]]+v2[m1,K[1]])), K[1], 1, m2])) (C[2])[m1]],m2,m1])//FullSimplify

p1[m1_,m2_]=(E^(InactiveIntegrate[-((Derivative[1,0][t1][K[1],m2])/(t1[K[1],m2]-v1[K[1],m2])), K[1], 1, m1]))

D[p1[m1,m2],m2](v2[m1,m2]-1+t1[m1,m2])-p1[m1,m2]D[1-t1[m1,m2],m2]//FullSimplify

(InactiveIntegrate[((Derivative[0,1][t1][s,m2] (Derivative[1,0][t1][s,m2]+Derivative[1,0][v2][s,m2])+(-1+t1[s,m2]+v2[s,m2])Derivative[1,1][t1][s,m2])/(((-1+t1[s,m2]+v2[s,m2]))^2)), s, 1, m1])//FullSimplify

Integrate[((Derivative[0,1][t1][s,m2] (Derivative[1,0][t1][s,m2]+Derivative[1,0][v2][s,m2])+(-1+t1[s,m2]+v2[s,m2])Derivative[1,1][t1][s,m2])/(((-1+t1[s,m2]+v2[s,m2]))^2)),{s,1,m1}]

D[Derivative[0,1][t1][s,m2]/(-1+t1[s,m2]+v2[s,m2]),s]

The above cancels so that the function only depends on m2.

The following line crashed...restate in simpler form...just an ODE in u[m1]=(1/(-1+t1[m1,m2]+v2[m1,m2]))v2 (1,0) derivative only so solve as v2[m1]

DSolve[compcond== 0, v2, {m1, m2}]

DSolve[(((Derivative[0,1][t1][m1,m2]-Derivative[0,1][v1][m1,m2]) Derivative[1,0][t1][m1,m2])/(((t1[m1,m2]-v1[m1,m2]))^2))-((Derivative[0,1][t1][m1,m2] (Derivative[1,0][t1][m1,m2]+v2'[m1]))/(((-1+t1[m1,m2]+v2[m1]))^2))-((Derivative[1,1][t1][m1,m2])/(t1[m1,m2]-v1[m1,m2]))+((Derivative[1,1][t1][m1,m2])/(-1+t1[m1,m2]+v2[m1]))== 0, v2, {m1}]

DSolve[(((Derivative[0,1][t1][m1,m2]-Derivative[0,1][v1][m1,m2]) Derivative[1,0][t1][m1,m2])/(((t1[m1,m2]-v1[m1,m2]))^2))-((Derivative[0,1][t1][m1,m2] (Derivative[1,0][t1][m1,m2]+v2'[m1]))/(((-1+t1[m1,m2]+v2[m1]))^2))-((Derivative[1,1][t1][m1,m2])/(t1[m1,m2]-v1[m1,m2]))+((Derivative[1,1][t1][m1,m2])/(-1+t1[m1,m2]+v2[m1]))==0,v2,{m1}]

DSolve[(((Derivative[0,1][t1][m1,m2]-Derivative[0,1][v1][m1,m2]) Derivative[1,0][t1][m1,m2])/(((t1[m1,m2]-v1[m1,m2]))^2))+Derivative[0,1][t1][m1,m2] u'[m1]-((Derivative[1,1][t1][m1,m2])/(t1[m1,m2]-v1[m1,m2]))+Derivative[1,1][t1][m1,m2]u[m1]==0,u,{m1}]

Solve[(1/(-1+t1[m1,m2]+v2))==((C[1])/(Derivative[0,1][t1][m1,m2]))+((InactiveIntegrate[((-Derivative[0,1][t1][K[1],m2] Derivative[1,0][t1][K[1],m2]+Derivative[0,1][v1][K[1],m2] Derivative[1,0][t1][K[1],m2]+t1[K[1],m2] Derivative[1,1][t1][K[1],m2]-v1[K[1],m2] Derivative[1,1][t1][K[1],m2])/(((t1[K[1],m2]-v1[K[1],m2]))^2)), K[1], 1, m1])/(Derivative[0,1][t1][m1,m2])),v2]//FullSimplify

The above is the necessary relationship between v2 and v1, note c1 is a function of m2.

Suppose v1=f1(t1) and v2=f2(t1) for known functions f1 and f2 (recall t2=1-t1 so it should be removed everywhere), then we solve...Replace t2 by 1-t1

DSolve[D[p[m1,m2],m1](f1[t1[m1,m2]]-t1[m1,m2])-p[m1,m2]D[t1[m1,m2],m1]==0,{p[m1,m2]},m1]

DSolve[D[p[m1,m2],m2](f2[t1[m1,m2]]-1+t1[m1,m2])-p[m1,m2]D[1-t1[m1,m2],m2]==0,{p[m1,m2]},m2]

compcond=(D[Log[(E^(InactiveIntegrate[((Derivative[1,0][t1][K[1],m2])/(f1[t1[K[1],m2]]-t1[K[1],m2])), K[1], 1, m1])) (C[1])[m2]],m1,m2]-D[Log[(E^(InactiveIntegrate[-((Derivative[0,1][t1][m1,K[1]])/(-1+f2[t1[m1,K[1]]]+t1[m1,K[1]])), K[1], 1, m2])) (C[2])[m1]],m2,m1])//FullSimplify

This derivative compatibility quantity must be 0 in order for there to be any solution.

DSolve[compcond==0,t1,{m1,m2}]

DSolve[D[p[m1,m2],m1](v1[m1,m2]-t1[m1,m2])-p[m1,m2]D[t1[m1,m2],m1]==0,{p[m1,m2]},m1]

DSolve[D[p[m1,m2],m2](v2[m1,m2]-1+t1[m1,m2])-p[m1,m2]D[1-t1[m1,m2],m2]==0,{p[m1,m2]},m2]

compcond=(D[Log[(E^(InactiveIntegrate[-((Derivative[1,0][t1][K[1],m2])/(t1[K[1],m2]-v1[K[1],m2])), K[1], 1, m1])) (C[1])[m2]],m1,m2]-D[Log[(E^(InactiveIntegrate[-((Derivative[0,1][t1][m1,K[1]])/(-1+t1[m1,K[1]]+v2[m1,K[1]])), K[1], 1, m2])) (C[2])[m1]],m2,m1])//FullSimplify

Set u=Log[p]

DSolve[{D[u[m1,m2],m1]==D[t1[m1,m2],m1]/(f1[t1[m1,m2]]-t1[m1,m2]),D[u[m1,m2],m2]==D[1-t1[m1,m2],m2]/(f2[t1[m1,m2]]-1+t1[m1,m2])},u,{m1,m2}]

DSolve[{D[u[m1,m2],m1]==r1[m1,m2],D[u[m1,m2],m2]==r2[m1,m2]},u,{m1,m2}]

DSolve[D[r1[m1,m2],m2]==D[r2[m1,m2],m1],r2,{m1,m2}]

D[((Derivative[1,0][t1][m1,m2])/(f1[t1[m1,m2]]-t1[m1,m2])),m2]==D[-((Derivative[0,1][t1][m1,m2])/(-1+f2[t1[m1,m2]]+t1[m1,m2])),m1]//FullSimplify

DSolve[((-(-1+f1'[t1]) t12 t11+(f1[t1]-t1) t112)/(((f1[t1]-t1))^2))==(((1+f2'[t1]) t12 t11-(-1+f2[t1]+t1) t112)/(((-1+f2[t1]+t1))^2)),f2,t1]

The above is the compatibility condition for f2 in terms of f1 as a function of t1. Here t12= Derivative[0,1][t1][m1,m2], t11= Derivative[1,0][t1][m1,m2], t112= Derivative[1,1][t1][m1,m2]. These are assumed to be known as is f1[t1], then f2[t2] is forced to be as above.

Next we do the case that v1 is a function only of m2 and v2 only of m1. The compatibility condition gives v2[m1] but there is dependency of m1.

DSolve[D[p[m1,m2],m1](v1[m2]-t1[m1,m2])-p[m1,m2]D[t1[m1,m2],m1]==0,{p[m1,m2]},m1]

DSolve[D[p[m1,m2],m2](v2[m1]-1+t1[m1,m2])-p[m1,m2]D[1-t1[m1,m2],m2]==0,{p[m1,m2]},m2]

compcond2=(D[Log[(((C[1])[m2])/(t1[m1,m2]-v1[m2]))],m1,m2]-D[Log[(((C[2])[m1])/(1-t1[m1,m2]-v2[m1]))],m2,m1])//FullSimplify

DSolve[compcond2==0,v2,m1]