If you end your printf statements with semicolons, then those characters will not appear.

The return value of a printf statement is NULL. When two printfs are juxtaposed in 2D input, the juxtaposition is interpreted as multiplication. NULL*NULL = NULL^2, which prints as ()^2. If there is an intervening semicolon, then they are not juxtaposed.

@John Fredsted It is not ad hoc: evalc assumes that variables are real unless otherwise specified.

In your first code, you were accidentaly assuming that the variables were real, but Maple was not. Maple could not put any variables on the diagonal in that case because they were not known to be real.

There is a three-dimensional figure called a Gomboc which was discovered with the help of Maple. See this Wikipedia page.

I'd recommend

1. The Art of Computer Programming, Volume 2: Seminumerical Algorithms, by Donald Knuth. This book has been translated into Spanish and many other languages.

2. Modern Computer Algebra, by Joachim von zur Gathen and Jurgen Gerhard.

3. Matrix Computations, by Gene H. Golub and Charles F. Van Loan.

You have 8 data points and 15 parameters. The number of data points must be greater than or equal to the number of parameters.

This is considered linear data fitting rather than nonlinear because the model function is linear in the parameters. The fact that it is nonlinear in the independent variables is irrelevant.

If you had enough data, the commands would be

x1_values:= < 0.1, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80 >:
x2_values:= < 1, 2, 3, 4, 5, 6, 7, 8 >:
x3_values:= < 11, 12, 13, 14, 15, 16, 17, 18 >:
x4_values:= < 10, 20, 30, 40, 50, 60, 70, 80 >:
y_values:= < 30, 40, 60, 70, 90, 120, 150, 200 >:
model:=
     a + b*x1 + c*x2 + d*x3 + e*x4 + f*x1^2 + g*x2^2 + h*x3^2 + i*x4^2 + j*x1*x2 +
     k*x1*x3 + l*x1*x4 + m*x2*x3 + n*x2*x4 + p*x3*x4
:

Statistics:-LinearFit(
     model,
     < x1_values | x2_values | x3_values | x4_values >,
     y_values,
     [x||(1..4)]
);

Note that the old stats package is deprecated and has been replaced by Statistics.

Your main problem is that you have used = for assignment statements when you should be using :=. Here's my Collatz procedure:

Collatz:= proc(val::posint)
local res:= val, k;
     for k while res <> 1 do
          if irem(res, 2, 'res') = 1 then res:= 6*res+4 end if;
     end do;
     k-1
end proc;

This uses the fact that irem also computes the quotient and returns it in the third argument, if present. You may ask "Why 6*res+4 instead of 3*res+1?" The reason is that I need to compensate for the fact that the integer quotient of the odd number has already been taken. So it's 3*(2*n+1) + 1, which is 6*n+4. The return value is the number of iterations required to get to 1.

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