Here is the code that lets you change the upper and lower bounds for the parameters and the numbers of increments for the paramters.

restart:
eq1:= diff(f(eta),eta$3) - a*diff(f(eta),eta$1)^2 + b*f(eta)*diff(f(eta),eta$2) = 0:
bc:= f(0)=0, D(f)(0) = 1+c*(D@@2)(f)(0), D(f)(8)=d:
sys:= unapply({eq1,bc}, a, b, c, d):
Sol:= (a,b,c,d)-> dsolve(sys(a,b,c,d), numeric):
Ranges:= [0..1, 0..1, 0..1, 0..1]: #Ranges for a, b, c and d, respectively
nIncrs:= [2,2,2,2]:  #Numbers of increments for a, b, c, d
Pts:= [seq([seq(Ranges[k], -`-`(op(Ranges[k]))/nIncrs[k])] , k= 1..nops(nIncrs))]:
It:= combinat:-cartprod(Pts):
A:= Matrix(`*`(nops~(Pts)[]), nops(Pts)+1):
for row while not It[finished] do
     V:= It[nextvalue]()[];
     A[row,..]:= < V, eval(diff(f(eta),eta$2), Sol(V)(0)) >
end do:

interface(rtablesize= infinity):
A;

@rashmi Any plots (of the same number of dimensions) can be combined with  plots:-display. For example, using the plot my Answer above,

P1:= plot([Re,Im](eval(rhs(Sol), {C1=1, C2=1})), tau= 0..2, axes= boxed):
P2:= plots:-textplot([1,20,"My solution"], font= [HELVETIC,BOLDITALIC,16]):
plots:-display([P1,P2]);

This example also shows how to set the text size in a textplot, the 16 being the size in the example.

You mentioned odeplots, so I thought that I should mention that odeplot is only for numeric dsolve solutions. The ODE solutions discussed above are symbolic, not numeric.

Instead of display, use plots:-display.

You could use solve with floating-point coefficients. But why do you want another way?

The problem with using combinat:-permute(1,1,2,2,3,3,4,4,5,5,6,6) is that it doesn't give an equiprobable space. To get an equiprobable space, use combinat:-cartprod([[$1..6] $ 2]).

Counts:= table(sparse):
It:= combinat:-cartprod([[$1..6]$2]):
while not It[finished] do
     V:= `+`(It[nextvalue]()[]);
     Counts[V]:= Counts[V]+1
end do:
eval(Counts);

Another way to get the distribution is with Statistics.

restart:
macro(St= Statistics):
DicePair:= `+`('St:-RandomVariable(ProbabilityTable([1/6 $ 6]))' $ 2):
St:-CDF(DicePair, 9) - St:-CDF(DicePair, 6);

The solution that you got is real for appropriate values of the constants. But it would be difficult to figure out those values. The trick is to use your own constants from the start:

ode:= diff(eta(tau), tau, tau)+(8/(4*tau^2+1)-32/(4*tau^2+1)^2)*eta(tau) = 0:
Sol:= dsolve({ode, eta(0)=C1, D(eta)(0)=C2}):

Now if you set C1 and C2 to real values, you should get a real solution. The problem now is that during numeric evaluation for plotting, there will be small spurious imaginary parts. These can be discarded by using Re. The following plot shows the plot of the solution and also shows that the imaginary part is effectively zero. Once you are confident that it is effectively zero, you can remove the Im from the plot.

plot([Re,Im](eval(rhs(Sol), {C1=1, C2=1})), tau= 0..2, axes= boxed);

What you want is signum(x). The sign command is just for polynomials, so sign(x) for unevaluated x is already 1 before x has been assigned a value.

restart:
eq1:= diff(f(eta),eta$3) - a*diff(f(eta),eta$1)^2 + b*f(eta)*diff(f(eta),eta$2) = 0:
bc:= f(0)=0, D(f)(0) = 1+c*(D@@2)(f)(0), D(f)(8)=d:
sys:= unapply({eq1,bc}, a, b, c, d):
Sol:= (a,b,c,d)-> dsolve(sys(a,b,c,d), numeric):
Incrs:= [.01, .1, 0.1, 0.1]: #Increments in a, b, c and d, respectively
nIncrs:= [2,2,2,2]:  #Numbers of increments
It:= combinat:-cartprod([seq([$1..nIncrs[k]]*Incrs[k], k= 1..nops(Incrs))]):
A:= Matrix(`*`(nIncrs[]), nops(Incrs)+1):
for row while not It[finished] do
     V:= It[nextvalue]();
     A[row,..]:= < V[], eval(diff(f(eta),eta$2), Sol(V[])(0)) >
end do:

This is an application of Huygens's gambler's ruin theorem. The probability of winning any one iteration is p = 5/12. I'll assume that you don't need any help computing that. We'll say that the player's starting "bank" is n[1] = 2 chips (two $50 bets), and, effectively, the casino's starting bank is also n[2] = 2 chips, because the game ends when either side has lost its starting bank. Then Huygens's formula for the probability that the casino loses its bank first is

(1 - (1/p-1)^n[1])/(1-(1/p-1)^(n[1]+n[2]));
                                 

eval(%, [p= 5/12, n[1]= 2, n[2]= 2]);
                              25/74

See this Wikipedia article

See the section entitled "Unfair coin flipping".

The inner loop is performed two times regardless of whether the n is actually accessed in the loop. Indeed, your double loop could just as well be expressed as

QQ:=Matrix([[3],[4],[1]]);

for m from 1 to 2 do
     to 2 do
          QQ:= (QQ*m)+QQ
     end do
end do:

Yes, it is easy. But is there any significance to the brackets? If there is no significance, then do

`and`(ListTools:-Flatten(L)[])

where L is your original list.

Brian,

My program can be easily modified to work with lists instead of sets, and real numbers instead of integers. Here is is:

S:= [3, 4, 5, 6, 8, 9, 28, 30, 35]:
SL:= [A,B,C,D,E,F,G,H,I]:
assign(Labels ~ (S) =~ SL); #Create remember table.
AllP:= combinat:-setpartition(S,3):
lnp:= evalf(ln((`*`(S[]))^(1/3))):

Var:= (P::({list,set}({list,set})))->
     evalf(`+`(map(b-> abs(ln(`*`(b[]))-lnp), P)[]))
:

Min:= proc(S::{list,set}, P::procedure)
local M:= infinity, X:= (), x, v;
     for x in S do
          v:= P(x);
          if v < M then  M:= v;  X:= x  end if
     end do;
     X
end proc:

ans:= Min(AllP, Var);
              
subsindets(ans, realcons, Labels);

Instead of using the linear model y = a + b*x, use y = b*x. If you're using Statistics:-LinearFit, then the command is

Statistics:-LinearFit(b*x, X, Y, x);

where X is your list or array of x data and Y is your list or array of y data. If you're using CurveFitting:-LeastSquares, then the command is

CurveFitting:-LeastSquares(X, Y, x, curve= b*x);

Surely this problem was contrived to point out the flaws in a floating-point arithmetic system. In particular, this problem looks contrived specifically to dupe Maple because the error occurs right at 10 decimal places, which is Maple's default. So I can't believe that it is truly a BIG problem for you.

The problem occurs just in the computation of

.5 - 1234567892;

This value requires 11 digits to represent, obviously, but you only have 10 digits to work with (by default). It is easy to increase the default:

Digits:= 11;

Now you will get accurate results for this computation.

Another option is to use exact arithmetic:

convert(.1234567891, rational, exact)*10^10 +1/2 - 1234567892;

Exact arithmetic works accurately regardless of the value of Digits.

I managed to make it into a simulation for any number of doors and still simplify and speed up the code some.

# Proc to make a uniform random selection from a set.
Rand:= (S::set)-> `if`(nops(S) < 2, S[], S[(rand() mod nops(S)) + 1]):

# Simulation for the n-door Monty Hall problem.
# Returns the number of trials the player wins.
MontyHall:= proc(
     {trials::posint:= 10000},
     {doors::posint:= 3},
     {switches::list(truefalse):= [true]},
     $
)
local
     k, wins:= 0, Doors,
     PrizeDoor, #The door with the prize behind it.
     PlayersDoor,  # The door the player picks.
     AllDoors:= {$1..doors}   
;
     to trials do
          Doors:= AllDoors;
          PrizeDoor:= Rand(Doors);
          PlayersDoor:= Rand(Doors);
          for k to doors-2 do
               # Monty opens a door, which is thus removed from the set
               # of available doors.
               Doors:= Doors minus {Rand(Doors minus {PrizeDoor,PlayersDoor})};       
               if switches[k] then
                    PlayersDoor:= Rand(Doors minus {PlayersDoor})
               end if
          end do;
          if PlayersDoor = PrizeDoor then  wins:= wins+1  end if
     end do;
     wins
end proc:     

MontyHall(trials= 100000, doors= 5, switches= [false,false,true]);
                             79904