The entropy function

entropy[x_,y_] := 
-(Which[x==0,0,True,x*Log[x]]+
Which[y==0,0,True,y*Log[y]]
);

maxentropy = 2/E;
maxplot =
Graphics[{Blue,Point[{1/E,1/E}]}]; 

cplots = 
ContourPlot[
{entropy[x,y]==(1/5)maxentropy,
entropy[x,y]==(2/5)maxentropy,
entropy[x,y]==(3/5)maxentropy,
entropy[x,y]==(4/5)maxentropy
},
{x,0,1},{y,0,1},
ContourStyle->{Blue}
];

Show[
cplots,maxplot
]

The entropy function with the probability constraint

probplot =
Plot[
1-x,
{x,0,1},
PlotStyle->Red
];

tangentpplot =
ContourPlot[
entropy[x,y]==entropy[1/2,1/2],
{x,0,1},{y,0,1},
ContourStyle->{Blue}
];

gpvector = 
Graphics[{Red,Arrow[{{1,1}/2,{1,1}/2+{1,1}/10}]}];
gevector =
Graphics[{Blue,Arrow[{{1,1}/2,{1,1}/2+{-1+Log[2],-1+Log[2]}/10}]}];


Show[
cplots,maxplot,tangentpplot,probplot,gpvector,gevector
]

The entropy function with the probability & average E constraints

Eave = -0.3;
Print["For (OverBar[E]) = ",Eave," :"];

{βsol} = NSolve[Eave == (-E^b + E^(-b))/(E^b + E^(-b)),b];
β = b/.βsol;
Print["    β = ",β];
Z = (E^β + E^(-β));
λ = 1 - Log[Z];
Print["    λ = ",λ];

jointplot = 
Plot[
(λ-β*Eave)/(λ-β) - x*(λ+β)/(λ-β),
{x,0,1},
PlotStyle->Red
];

p0 = E^β/Z;
p1 = E^(-β)/Z;

tangentjplot =
ContourPlot[
entropy[x,y]==entropy[p0,p1],
{x,0,1},{y,0,1},
ContourStyle->{Blue}
];

gjvector = 
Graphics[
{Red,
Arrow[{{p0,p1},{p0,p1}+{λ+β,λ-β}/10}]
}
];
gevector = 
Graphics[
{Blue,
Arrow[{{p0,p1},{p0,p1}+{-1-Log[p0],-1-Log[p1]}/10}]
}
];


Show[
cplots,maxplot,tangentpplot,tangentjplot,jointplot,gjvector,gevector
]