@Carl Love 

Dear Carl

I try to plot 

f:=epsilon=(lambda1+sqrt(-lambda2^2+Y*(lambda3^2+(lambda4-lambda5*Y)^2)))^2:

as 3d

implicitplot3d(f, Y = 0 .. 0.4,epsilon=0..20,Delta=-25*Pi..25*Pi,labels=[Y,E1,Delta],tickmarks=[3,3,3]);
can u plz check if it is correct

and the contour

with(Student[Calculus1]):contourplot(f,Delta =-25*Pi .. 25*Pi,Y=0..1,contours=[10],axes=boxed,thickness=2,color=black,font=[1,1,18],tickmarks=[5, 5],linestyle=1,view=[0..1,0..0.6]);

@Carl Love 

Thanks Carl

Deos it make a difference if I plot Y agianst (epsilon^2/(4Pi^2))

@Christian Wolinski 

many thankx for ur reply...

why you are calculating the roots..

Im really grateful...

I excpect optical bistability shape...

but many thanks

thanks for ur code

------------------------- Defining the nature of the variables used ----------------------

epsilon:=epsilon1/(2*Pi):gamma1:=8*Pi:gamma2:=0.002*Pi:x1:=100*Pi:omega2:=200*Pi:delta1:=0.2*Pi:Delta:=25*Pi:G:=20*Pi:Omega:=0:

lambda1:=G*Omega*gamma1/(2*(0.25*gamma1^2+Delta^2)):lambda2:=G*Omega*Delta/(0.25*gamma1^2+Delta^2):lambda3:=gamma1+lambda1:lambda4:=delta1-G^2*Delta/(0.25*gamma1^2+Delta^2):lambda5:=2*x1^2*omega2/(omega2^2+gamma2^2):c:=a^2+b^2:
f:=epsilon-(lambda1-I*lambda2)=(a+I*b)(lambda3+I*(lambda4+lambda5*(c))):
f:=convert(epsilon-(lambda1-I*lambda2)=(a+I*b)*(lambda3+I*(lambda4+lambda5*(a^2+b^2))),rational):
f1:=(lhs-rhs)(evalc(f)):

B:=coeff(f1,I):    

A:=expand(f1-B*I):
    
F:=t->fsolve({eval(A,epsilon=t),B}):
    
plot('eval(a^2+b^2,F(t))', t=0..5, labels=[epsilon,a^2+b^2], view=[0..5,0..0.5], labelfont=[times,16],color=black,axes=boxed,font=[1,1,18],tickmarks=[4,2],thickness=2);
Warning,  computation interrupted

@acer 

thanks.. but how can i know the  approximation function lambda2=ao+a1lambda1+a2lambda1^2

@acer 

thanks for ur answer,,

what i mean is for example we have y=x^2 the approximation function is : y=ao+a1x+a2x^2+a3x^3 

all we need 4 pair (x,y) which are (0,0),(1,1),(2,4),(3,9)

and then subsitue the points in the approximation function to find ao,a1,a2,a3

in the plot below all I need 4 point from the figure and substitue in y=ao+a1x+a2x^2+a3x^3

contourplot(4*lambda2*result^2/(Pi*(lambda2+1)^2)-lambda1,lambda1=0..1,lambda2=0..1,contours=[0],axes=boxed,title=tit,titlefont=[SYMBOL,16],thickness=1,color=black,font=[1,1,18],tickmarks=[2, 4],linestyle=1,view=[0.002..1,0.002..1]);

@Ramakrishnan 

how did u got it

mamy thanks

@acer 

k1:=sqrt(1-lambda1^2):
result:=int(sqrt(1-k1^2*sin(x)*sin(x)),x=0..Pi/2):

@Rouben Rostamian  

Many thanks that is useful 

@Rouben Rostamian  

In the Deplot3D it must have initial conditions, but the example I have is plot of the phase portraits for the system where the inc are unkown

                 /          (1/2)\         (1/2)         
                 |/    2   \     | /  1   \      / 2    \
      2 EllipticE||- ------|     | |------|      \a  - 1/
                 ||   2    |     | | 2    |              
                 \\  a  - 1/     / \a  - 1/              
                      /          (1/2)\         (1/2)
                      |/    2   \     | / 2    \     
           2 EllipticE||- ------|     | \a  - 1/     
                      ||   2    |     |              
                      \\  a  - 1/     /              
                  /                 (1/2)\         (1/2)
                  | (1/2) /    1   \     | / 2    \     
       2 EllipticE|2      |- ------|     | \a  - 1/     
                  |       |   2    |     |              
                  \       \  a  - 1/     /              

@ThU 

that is easy many thanks

I was thinking of (wich is work)

eq1:=u=(3*x-1)*(1-x):
eq2:=v=2*x*(1-2*x):
U:=[solve(eval(eq1),x)]: 
vua:=eval(solve(eq2,v)); 
V:=eval~(vua,x=~U): 
P3:=plot(V,u=-1..1,view=[-1..1,-2..2],color=red,labels=[u,v],axes=boxed,numpoints=9000,linestyle=1,font=[1,1,18],tickmarks=[3,4],thickness=2);

@Carl Love 

thanks .. true the roots (for these fixed parameters) are between 0..1 & 8..10.

and in general (changing the parameters) the roots will be different.

many thanks indeed. 

@tomleslie 

1- the plots of the +ve roots aganst arbitrary so(4..5) is done by carl.

2 - want to subsitute these positive roots in the expressions xs,zs also for arbitrary so for example: 

if so=4.3 then we sebstitue it in Q --> we get two +ve roots and  then subsitute them in xs hence, zs ..etc then plot xs against so, plot zs against so.

many thanks for ur patients.