An example would help tremendously.

An example would help tremendously.

@asfn In y appear 3 square roots: sqrt(12.25 - b^2), sqrt(2.1025-b^2), and sqrt(1-b^2). This reveals where changes are likely.

Without attempting a proof that my statement about y being purely imaginary between 1 and sqrt(2.1025) you can test the validity by evaluating y at points chosen randomly:

y:=subs(`βb`=b,y);
rnd:=RandomTools:-MersenneTwister:-GenerateFloat;
#Testing rnd
rnd();

#Now evaluating y at 10 points chosen randomly in the strip mentioned:

to 10 do evalf(eval(y,{b=1+rnd()*(sqrt(2.1025)-1),lambda=1+rnd()})) end do;

Similarly, you can test the claim that y is real for sqrt(2.1025) < b < sqrt(12.1025) by

to 10 do evalf(eval(y,{b=sqrt(2.1025)+rnd()*(sqrt(12.25)-sqrt(2.1025)),lambda=1+rnd()})) end do;

Now the disturbing question is: What did implicitplot actually do? It should have painted the whole strip 1 < b < sqrt(2.1025) red, but produced some nice looking curves!

@asfn In y appear 3 square roots: sqrt(12.25 - b^2), sqrt(2.1025-b^2), and sqrt(1-b^2). This reveals where changes are likely.

Without attempting a proof that my statement about y being purely imaginary between 1 and sqrt(2.1025) you can test the validity by evaluating y at points chosen randomly:

y:=subs(`&beta;b`=b,y);
rnd:=RandomTools:-MersenneTwister:-GenerateFloat;
#Testing rnd
rnd();

#Now evaluating y at 10 points chosen randomly in the strip mentioned:

to 10 do evalf(eval(y,{b=1+rnd()*(sqrt(2.1025)-1),lambda=1+rnd()})) end do;

Similarly, you can test the claim that y is real for sqrt(2.1025) < b < sqrt(12.1025) by

to 10 do evalf(eval(y,{b=sqrt(2.1025)+rnd()*(sqrt(12.25)-sqrt(2.1025)),lambda=1+rnd()})) end do;

Now the disturbing question is: What did implicitplot actually do? It should have painted the whole strip 1 < b < sqrt(2.1025) red, but produced some nice looking curves!

G1:=Matrix(101,93,(i,j)->evalf(sin(i/7)+cos(j/5)),datatype=float[8]):

plots:-surfdata(convert(G1,listlist));

plots:-surfdata(Array(G1));

G1:=Matrix(101,93,(i,j)->evalf(sin(i/7)+cos(j/5)),datatype=float[8]):

plots:-surfdata(convert(G1,listlist));

plots:-surfdata(Array(G1));

@zhkhan There appears a G, which apparently is just a name.

@zhkhan There appears a G, which apparently is just a name.

@zhkhan I'm not sure what you mean.

But the two larger roots start their existence at the following value for l: 

fsolve({diff((rhs-lhs)(A),sigma)=0,(rhs-lhs)(A)=0},{sigma,l=3.5});

 Output: {l = 3.531408488, sigma = -32.34143629}

I just used the default setting of Digits = 10.

Both roots can be found by RootFinding:-NextZero:

r1:='r1': r2:='r2':
for l from 3.7 by .1 to 5 do
   r1:=RootFinding:-NextZero(unapply((rhs-lhs)(A),sigma),-120);
   if r1<>FAIL then r2:=RootFinding:-NextZero(unapply((rhs-lhs)(A),sigma),r1) end if;
   print(l,r1,r2);
end do:
l:='l':

There will be numerical problems when the first of the two roots is close to the singularity of rhs(A), which happens for larger values of l.

@zhkhan I'm not sure what you mean.

But the two larger roots start their existence at the following value for l: 

fsolve({diff((rhs-lhs)(A),sigma)=0,(rhs-lhs)(A)=0},{sigma,l=3.5});

 Output: {l = 3.531408488, sigma = -32.34143629}

I just used the default setting of Digits = 10.

Both roots can be found by RootFinding:-NextZero:

r1:='r1': r2:='r2':
for l from 3.7 by .1 to 5 do
   r1:=RootFinding:-NextZero(unapply((rhs-lhs)(A),sigma),-120);
   if r1<>FAIL then r2:=RootFinding:-NextZero(unapply((rhs-lhs)(A),sigma),r1) end if;
   print(l,r1,r2);
end do:
l:='l':

There will be numerical problems when the first of the two roots is close to the singularity of rhs(A), which happens for larger values of l.

@zhkhan You may want to take a look at ExportMatrix.

@zhkhan You may want to take a look at ExportMatrix.

Since you are solving the ode numerically it may be a better idea to keep the second order equation (i.e. Newton's 2nd law). Your first order equation is obtained by multiplying the 2nd order equation by diff(theta(t),t) and integrating. This equation is not equivalent to the original.

Since you are solving the ode numerically it may be a better idea to keep the second order equation (i.e. Newton's 2nd law). Your first order equation is obtained by multiplying the 2nd order equation by diff(theta(t),t) and integrating. This equation is not equivalent to the original.

Could you provide the actual code?