# Some simple Calculations of the bifurcation diagram for the Lorenz system
> restart;
> with(linalg):
> with(plots):
> vals:= sigma=10,b=8/3;
> allvals:=vals,r=28;
> xdot := sigma*(y-x);
> ydot:= r*x-y-x*z;
> zdot:= x*y-b*z;
> equil:=solve({xdot,ydot,zdot},{x,y,z});
> eq1:=equil[1];
> eq2:=allvalues(equil[2])[1];
> Jac:= jacobian([xdot,ydot,zdot],[x,y,z]);
> J_0:= subs(eq1,evalm(Jac));
> l:= eigenvalues(J_0);
> ll:= subs(vals,{l});
> plot1:= plot({ll[1],ll[2],ll[3]},r=0..20, color=[black,blue,orange]):
> display(plot1);
> J1:= subs(eq2,evalm(Jac));
> p2:= collect(charpoly(J1,lambda), lambda);
> p2:= subs(vals,%);
> rts := solve(p2,lambda):
> plot2:= plot({Re(rts[1]),Re(rts[2]),Re(rts[3])},r=1..30, color = [red,yellow,green],thickness=2):
> display(plot1,plot2,view=[23..30,-1..1]);
> plot3:= plot({Im(rts[1]),Im(rts[2]),Im(rts[3])},r=1..10, color = [red,yellow,green],thickness=2):
> display(plot3,view=[0..5,-5..5]);