This is not "solving" the conjecture; it is verifying it.

The following could be done faster by sieving for the primes, at the expense of a few more lines of code.

Goldbach:= proc(n::even)
local ct:= 0, p:= 2;
     while p <= n/2 do
          if isprime(n-p) then ct:= ct+1 end if;
          p:= nextprime(p)
     end do;
     ct
end proc:
    
Goldbach(10^7);
                   38807

X:= proc(n)
option remember;
     n*(thisproc(n-3)+thisproc(n-2)+thisproc(n-1))
end proc:
(X(0),X(1),X(2)):= (0,1,2):

length(X(2013));
                              5782

There may be some number-theoretic subtlety that can be applied to this problem. However 2013^1102 is not a big number by Maple standards, and this problem can be done in a few seconds by straightforward computation.

It is not clear to me that every number is attracted to one of the three orbits that you gave. The following procedure only terminates if one of those orbits is reached.

f:= x-> `if`(x::even, x/2, 3*x-1):

F:= proc(X::posint)
     local x:= X;
     do
          if x in {1,5,17} then return x end if;
          x:= f(x)
     end do
end proc:

F(2013^1102);
                              1

When applied to functions (in the Maple sense of that word), combine does the opposite of expand. So, if you agree with the output of expand(cos(x+y)), then the logical extension is cos(2*x) = cos(x+x). So expand(cos(2*x)) is 2*cos(x)^2 - 1 = 1 - 2*sin(x)^2. The combine is just applying that in reverse.

It is a bug, not anything wrong with your code. It gets stuck in an infinite loop. I haven't figured out a workaround yet. If you simplify the model to f:= a+b*x^c, then it works immediately.

See ?Digits and ?Rounding . For chopping:

Digits:= 4:  Rounding:= 0:

For rounding:

Digits:= 4:  Rounding:= nearest:

Go to the Maple Applications Center to get the Orthogonal Expansions package. Click on "Download attached file". Follow the trivial installation instructions in the ReadMe file.

Now you can easily find Fourier series in Maple. (Note that Pi is spelled with a capital P in Maple.)

restart:
f:= x-> min(abs(x), Pi/2):
FS:= OrthogonalExpansions:-FourierSeries(f(x), x= -Pi..Pi, infinity);

plot(f(x), x= -Pi..Pi);

plot([seq(eval(FS, infinity= n),  n= [1,2,3,5,6,30])], x= -Pi..Pi);

You will gain understanding of the order in which things are done if you set printlevel:= 15 before running your loop. You are reusing k for every i and j. So only the last (i,j) combination is saved for each k. If you don't understand that, I can explain more. Let me know.

What you want is this:

for i to 3 do
     for j to 3 do
          k:= 3*(i-1)+j;
          a[k]:= (i,j)
     end do
end do;

eval(a);

Note that table entries are not necessarily stored in any order.

Your nested for-do loops can also be done by nested seq loops and a table constructor:

a:= table([seq(seq(3*(i-1)+j = (i,j), j= 1..3), i= 1..3)]);

Try

subs(sqrt(c)= omega*sqrt(m), xsol);
simplify(%);

Try this:

restart;
eta:= (x,y)-> y*sqrt(U*v/x);
alias(eta= eta(x,y));
Phi:= (x,y)-> sqrt(U*v*x)*f(eta);
e3:= vx = diff(Phi(x,y),x);
e4:= vy = diff(Phi(x,y),y);

sort([B||(1..3)], (A,B)-> (A[2]-A[1]) > (B[2]-B[1]))[1];

This looks at the width of the interval, not just the upper limit of the interval.

The sort is builtin and very fast. But if you were doing a very large number of intervals (millions probably), it may be faster to code a linear search for the max.

It's very easy:

L:= [[0,2,3,0,0], [2,3,0,0,0], [1,1,2,2,0]]:
add(1/mul(x!, x= LL), LL= L);

Also, it is not necessary that the lists be the same length. So there is no need to pad with zeros.

Clearly there is a bug here. Here is a workaround: Bring the w inside and switch the order of summation:

sum(sum(binomial(2*q,n-1)*z^q*w^n/q!, n=1..infinity),q=0..infinity);

Thanks for a quite interesting problem.

Unfortunately, dsolve doesn't allow explicit parameters for BVPs. But we can trick it into taking them by calling dsolve within a procedure whose parameters are the parameters that we want. I take the first three of your four bc and use them in dsolve. Then I make an implicitplot using the fourth bc.

Eq:= diff(f(y),y$3) + f(y)*diff(f(y),y$2) - diff(f(y),y$1)^2 -
     S*(diff(f(y),y$1) + y/2*diff(f(y),y$2))=0:
bc:= D(f)(0)=1, f(0)=0, (D@@2)(f)(b)=0:

F:= proc(_b, _S)
local Sol:= dsolve(eval({Eq,bc}, {S= _S, b= _b}), f(y), numeric);
     eval(f(y), Sol(_b)) - _S*_b/2
end proc:
     
plots:-implicitplot(F, 1..3, 1..3);

It will take some guess work adjusting the ranges to get the plot that you want. There are numerous options to implicitplot that can be used for refinement. A plot3d of F may help you find the level curve F=0. Note that b is the first argument of F, so b is on the horizontal axis.

My oh my, you've certainly made this far more complicated and redundant than it needs to be. You don't need odeplot, pointplot, or display; simply plot will do. Simplify your code to this:

restart:
A := 0.2e-1: B := 10^(-5): k := 0:
eq1 := diff(X(t), t) = -(A+B*X(t))*X(t):
ic[1] := X(365*k) = 1000:
s[1] := dsolve({eq1, ic[1]}, X(t), range = 0 .. .365, numeric);
pt[1] := eval([t, X(t)], s[1](0));
for k to 3 do
     tk := 365*k;
     A := rhs(s[k](tk)[2]);
     ic[k+1] := X(tk) = 500.*A;
     s[k+1] := dsolve({eq1, ic[k+1]}, X(t), range = tk .. 2*tk, numeric);
     pt[k+1] := eval([t, X(t)], s[k+1](tk));
end do;

plot([seq(pt[k], k= 1..4)], axes= boxed);