> 
# lorenz.ms 
# Phase Portraits For Lorenz Equations
# 
# Each call to lorenz can produce 5 plots. Comment out the print
# statements for the plots that you don't want.
# 
> lorenz := proc(sigma, b, r, ts2, te2, ts1, te1, steps, fnameprefix)
>    local de1, de2, de3, desys, g, k, s, OPT, caption, xt,
>       dim, tp, xp, yp, zp, tt, dt;
> 
>    de1 := diff(x(t),t) = sigma*(-x(t)+y(t));
>    de2 := diff(y(t),t) = r*x(t) - y(t) - x(t)*z(t);
>    de3 := diff(z(t),t) = -b*z(t) + x(t)*y(t);
>    desys := {de1,de2,de3,x(0)=1,y(0)=1,z(0)=1}, [x(t),y(t),z(t)];
>    
>    dim := steps;
>    tp := array(1..dim);
>    xp := array(1..dim);
>    yp := array(1..dim);
>    zp := array(1..dim);
> 
>    # Solve the system and save results in arrays for plotting
>    # only save long term behaviour.
>   
>    k := 0;
>    dt := (te2-ts2)/(steps-1);
>    g := dsolve(desys, type=numeric,'maxfun'=-1);
>    for tt from ts2 to te2 by dt do
>       s := g(tt);
>       k := k + 1;
>       tp[k] := tt;
>       xp[k] := rhs(s[2]); yp[k] := rhs(s[3]); zp[k] := rhs(s[4]);
>    od;
>    dim := k;
> 
>    # plot some projections of phase space trajectories
> 
>    caption := cat(`r = `,convert(r,string));
>    OPT := axes=BOX, title = caption;
>    #interface(plotoutput=fnameprefix.`xy.ps`);
>    print(plot( [seq( [xp[k],yp[k]], k=1..dim)],labels=[`x`,`y`],OPT));
>    #interface(plotoutput=fnameprefix.`xz.ps`);
>    print(plot( [seq( [xp[k],zp[k]], k=1..dim)],labels=[`x`,`z`],OPT));
>    #interface(plotoutput=fnameprefix.`yz.ps`);
>    print(plot( [seq( [yp[k],zp[k]], k=1..dim)],labels=[`y`,`z`],OPT));
>    #interface(plotoutput=fnameprefix.`xyz.ps`);
>    print(plots[spacecurve]([seq( [xp[k],yp[k],zp[k]], k=1..dim)],
>       labels=[`x`,`y`,`z`],shading=NONE,OPT));
> 
>    # resolve the system on the larger interval ts1<=t<=te1
>    # to graph x(t) against t
> 
>    g := dsolve(desys, type=numeric,'maxfun=-1');
>    #interface(plotoutput=fnameprefix.`tx.ps`); 
>    xt := plots[odeplot](g, [t,x(t)], ts1..te1,
>       numpoints=2000,title = cat(caption,` (x versus t)`));
>    print(plots[display]([xt],labels=[`t`,`x`],OPT));
> end:


# 
# TRY ONLY ONE CALL TO LORENZ AT A TIME: COMMENT OUT THE OTHERS
# 
> #interface(plotdevice=postscript);


# 
# Nonsymmetric projection in phase space
# 
# Simple periodic motion,  r = 225
> lorenz(10, 8/3, 225.0, 3.0, 3.5, 0.0, 5.0, 501,`L225`);


# 
# Symmetric projections in phase space
# 
# Simple periodic motion,  r = 475  
> # lorenz(10, 8/3, 475.0, 3.0, 3.34, 0.0, 5.0, 501, `L475`);  


# Period doubled periodic motion,  r = 163.5  
> # lorenz(10, 8/3, 163.5, 6.0, 7.2, 0.0, 10.0, 1001, `L1635`);  


# Another period doubling, r = 93  
> # lorenz(10, 8/3, 93.0, 10.0, 12.5, 0.0, 13.0, 1001, `L93`);  


# Another period doubling, r = 69.75  
> # lorenz(10, 8/3, 69.75, 12.0, 15.6, 5.0, 15.6, 1001, `L6975`);  


# Chaotic, strange attractor, r = 28  
> # lorenz(10, 8/3, 28.0, 10.0, 40.0, 0.0, 30.0, 3001, `L28`);


>