@Carl Love 

Many thanks Carl

if I want to use the +ve roots to calculate 

xs:=(so-u)*d*(a1+u)/(m1*u);

zs:=((m2-d2)*xs-a2*d2)*(a3+ys)/(m3*(a2+xs));

where u are the +ve roots of Q

I try to do

g1:=eval(subs(u=S[1],xs)); and so the other roots but how to put these expresions in loop cuz so(from 4..5)

REGARDS

@tomleslie 

Dear i try to make it clear as i understand see below (need ur comments)

Plot the roots S[1, 2, 3] against continous value of so,

Plot xs &zs against so
restart:
assume(so,real):
m1:=5:m2:=2:m3:=1.5:a1:=0.16:a2:=0.45:a3:=0.833:d1:=0.25:d2:=0.1:d3:=0.075:
ys:=a3*d3/(m3-d3):
d:=0.5:so:=4.9:
w3:=d*(m1-d1):w2:=d*d1*so-d*m1*so-2*a1*d*d1+a1*d*m1-a2*m1^2+a2*d1*m1+m1*m2*ys:w1:=2*a1*d*d1*so-a1*d*m1*so-a1^2*d*d1+a1*a2*d1*m1+a1*m1*m2*ys:wo:=a1^2*d*d1*so:
Q:=s->w3*s^3+w2*s^2+w1*s+wo:
for so from 4 to 5 do # s0 from 4 to 5(continous value)
sol:=evalf(solve(Q(s),s)): S:=array([],1..3): S[1]:=sol[1];S[2]:=sol[2];S[3]:=sol[3];

xs:=(so-u)*d*(a1+u)/(m1*u);zs:=((m2-d2)*xs-a2*d2)*(a3+ys)/(m3*(a2+xs));
g1:=eval(subs(u=S[1],xs));g2:=eval(subs(u=S[2],xs));g3:=eval(subs(u=S[3],xs));
k1:=eval(subs(u=S[1],zs));k2:=eval(subs(u=S[2],zs));k3:=eval(subs(u=S[3],zs));
end do;
with(plots):pointplot([[so,k1],[so,k2],[so,k3]],axes=boxed);# all value of [so,S[1,2,3] and [so,xs],[so,zs]

as u c the functions xs,zs depend on the three roots.

Any improvements of the above code (im sure its not right)

for fixed so=4

there are 2 +ve roots

restart:
assume(so,real):
m1:=5:m2:=2:m3:=1.5:a1:=0.16:a2:=0.45:a3:=0.833:d1:=0.25:d2:=0.1:d3:=0.075:
ys:=a3*d3/(m3-d3):
d:=0.5:so:=4:
w3:=d*(m1-d1):w2:=d*d1*so-d*m1*so-2*a1*d*d1+a1*d*m1-a2*m1^2+a2*d1*m1+m1*m2*ys:w1:=2*a1*d*d1*so-a1*d*m1*so-a1^2*d*d1+a1*a2*d1*m1+a1*m1*m2*ys:wo:=a1^2*d*d1*so:
Q:=s->w3*s^3+w2*s^2+w1*s+wo:Q(s):
sol:=evalf(solve(Q(s),s)): S:=array([],1..3): S[1]:=sol[1];S[2]:=sol[2];S[3]:=sol[3];

so how can we do it in for loop (so between 4 ..5) eact point calculate the roots then choose the +ve roots

finaly plot (+ve root,so)

@tomleslie 

for fixed so=4

there are 2 +ve roots

restart:
assume(so,real):
m1:=5:m2:=2:m3:=1.5:a1:=0.16:a2:=0.45:a3:=0.833:d1:=0.25:d2:=0.1:d3:=0.075:
ys:=a3*d3/(m3-d3):
d:=0.5:so:=4:
w3:=d*(m1-d1):w2:=d*d1*so-d*m1*so-2*a1*d*d1+a1*d*m1-a2*m1^2+a2*d1*m1+m1*m2*ys:w1:=2*a1*d*d1*so-a1*d*m1*so-a1^2*d*d1+a1*a2*d1*m1+a1*m1*m2*ys:wo:=a1^2*d*d1*so:
Q:=s->w3*s^3+w2*s^2+w1*s+wo:Q(s):
sol:=evalf(solve(Q(s),s)): S:=array([],1..3): S[1]:=sol[1];S[2]:=sol[2];S[3]:=sol[3];
0.008805705582
8.229389883
-0.07437287268

@tomleslie 

the aim is to plot only the +ve root against so

im thinking for do loop so from 4..5 and solve the cubic and choose the +ve roots so we get points (+ve root,so)

bucuse so must be between 4 to 5

best regards

@tomleslie 

Many thanks.

How about if we make table (+ve roots of Q(s),so) and plots the points

@Kitonum 

Thanks for ur clarifications

@Kitonum 

the warning is about the singularity,, how can we avoid if

@Carl Love 

Sorry If I do waste your time ...

many thanks

@Carl Love 

To calculate 

DeltaE:=2*(j*(spin*neighbours)+B*spin: p=exp(-DeltaE):

and 

calcultate the esponential 

p=exp(-DeltaE):

and decide which transition will occur

trans:=(rand(N)<trans)*(rand(N)<0.1)*-2+1:

perform the transition

spin:=spin*trans:

sum up our variables:

M:=sum(sum(spin)): E:=-sum(sum(DeltaE))/2:

display the current state of the system:

image((spin+1)*128);

@Carl Love @Doug Meade

I dont skip the line 

N:=3:B:=0:
j:=rand(0..1)()+10^(-10):
Spin := N -> (-1)^~(LinearAlgebra:-RandomMatrix( N, N, generator=rand(N) ));
spin := Spin(N):
N -> ~[^](-1, $,

LinearAlgebra:-RandomMatrix(N, N, generator = rand(N)))
for i from 1 to 2 do neighbours:=< spin[..,-1] | spin[..,1..-2] >+ < spin[..,2..-1] | spin[..,1] >+< spin[-1,..],spin[1..-2,..]>+ < spin[2..-1,..] , spin[1,..] >: DeltaE:=2*(j*(spin*neighbours)+B*spin: p=exp(-DeltaE): trans:=(rand(N)<trans)*(rand(N)<0.1)*-2+1: spin:=spin*trans: M:=sum(sum(spin)): E:=-sum(sum(DeltaE))/2: end do;

@Carl Love 

N:=3:B:=0:
j:=rand(0..1)()+10^(-10):
Spin := N -> (-1)^~(LinearAlgebra:-RandomMatrix( N, N, generator=rand(N) )); spin := Spin(N):
N -> ~[^](-1, $,

LinearAlgebra:-RandomMatrix(N, N, generator = rand(N)))
for i from 1 to 2 do neighbours:=< spin[..,-1] | spin[..,1..-2] >+ < spin[..,2..-1] | spin[..,1] >+< spin[-1,..],spin[1..-2,..]>+ < spin[2..-1,..] , spin[1,..] >: DeltaE:=2*(j*(spin*neighbours)+B*spin: p=exp(-DeltaE): trans:=(rand(N)<trans)*(rand(N)<0.1)*-2+1: spin:=spin*trans: M:=sum(sum(spin)): E:=-sum(sum(DeltaE))/2: end do;

Plz can u run the above and what is missing

@Doug Meade 

what in maple can I write this

image(C) displays the data in array C as an image. Each element of C specifies the color for 1 pixel of the image. The resulting image is an m-by-n grid of pixels where m is the number of columns and n is the number of rows in C. The row and column indices of the elements determine the centers of the corresponding pixels.

image((spin+1)*128);

@Doug Meade 

N:=3:B:=0:
j:=rand(0..1)()+10^(-10):
Spin := N -> (-1)^~(LinearAlgebra:-RandomMatrix( N, N, generator=rand(N) )); spin := Spin(N):
N -> ~[^](-1, $,

LinearAlgebra:-RandomMatrix(N, N, generator = rand(N)))
for i from 1 to 2 do neighbours:=< spin[..,-1] | spin[..,1..-2] >+ < spin[..,2..-1] | spin[..,1] >+< spin[-1,..],spin[1..-2,..]>+ < spin[2..-1,..] , spin[1,..] >: DeltaE:=2*(j*(spin*neighbours)+B*spin: p=exp(-DeltaE): trans:=(rand(N)<trans)*(rand(N)<0.1)*-2+1: spin:=spin*trans: M:=sum(sum(spin)): E:=-sum(sum(DeltaE))/2: end do;

@Doug Meade 

N:=1000:
Spin := N -> (-1)^~(LinearAlgebra:-RandomMatrix( N, N, generator=rand(N) ));
N -> ~[^](-1, $,

LinearAlgebra:-RandomMatrix(N, N, generator = rand(N)))
neighbours:=< spin[..,-1] | spin[..,1..-2] >+ < spin[..,2..-1] | spin[..,1] >+< spin[-1,..],spin[1..-2,..]>+ < spin[2..-1,..] , spin[1,..] >;

Error, (in rtable/Sum) invalid arguments