Everything on the left side of eq2 cancels, leaving just 0:

simplify(eq2);

So there is no C left to solve for.

Change your plot command from

plot([f(r)[1],f(r)[2],-18..18]);

to

plot([f(r)[1], f(r)[2], u=-18..18]);

Your code can be vastly simplified. I'll do that later.

Give the command

interface(rtablesize= 12);

Then you will be able to see the contents of your Matrix. I picked 12 because that's the largest dimension of your Matrix.

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);