(* Define the bootstrap matrix. *)
M={{1,xp[1,0],0},{xp[1,0],xp[2,0],I/2},{0,-I/2,1/3*(4*ene-xp[2,0])}}
(*xp[1,0]=0*)

(* Define the energy eigenvalue search space. *)
ene1=0;
ene2=1;
Nene=100;
energylist=Subdivide[ene1,ene2,Nene];

(* Choose the variable that you want to minimize. *)

minimizedvariable=xp[2,0]

(* Choose the solver and the tolerance. *)


method={"CSDP",Tolerance->10^(-10)}
method={"MOSEK",Tolerance->10^(-6)}
method={"DSDP",Tolerance->10^(-10)}


(* The main code for the bootstrap analysis. *)

dataMin={};
dataMax={};
Do[
bootstrapMatrix=M/.ene->i;
bootstrapvariables=Variables[bootstrapMatrix];

constraints= VectorGreaterEqual[{bootstrapMatrix,0}, {"SemidefiniteCone",Length[bootstrapMatrix[[1]]]}];

tmp=SemidefiniteOptimization[minimizedvariable,{constraints},bootstrapvariables,Method->method];

If[NumberQ[tmp[[1,2]]],AppendTo[dataMin,{i,tmp}]];

tmp=SemidefiniteOptimization[-minimizedvariable,{constraints},bootstrapvariables,Method->method];

If[NumberQ[tmp[[1,2]]],AppendTo[dataMax,{i,tmp}]];
,{i,energylist}]

(* Plot the allowed region: energy vs. x^2. *)

ListPlot[{Table[{i[[1]],minimizedvariable/.i[[2]]},{i,dataMax}],Table[{i[[1]],minimizedvariable/.i[[2]]},{i,dataMin}]}]