@Markiyan Hirnyk 

Using my original code, if I increase the initial point count using the 'number' option, as in

DirectSearch[GlobalSearch](abs(eqn1), assume=posint, number=500):

I get 249 solutions - I suspect there are many, many more :-(

@Rouben Rostamian  

As my original response said, isolve() provides an accurate solution in terms of eight parameters and I actually reached the same point that you did, ie

5 _Z1 + 10 _Z2 + 25 _Z3 + 35 _Z4 + 25 _Z5 + 40 _Z6 + 5 _Z7 + 10 _Z8 <666400

in which _Z1, _Z2, _Z3, _Z4, _Z5, _Z8  are required to be positive (not arbitrary, as you stated) integers:  _Z6, and  _Z7 may be negative, but lower-bounded by the conditions that f,g must be positive, ie

 _Z6 >  -(1/3)*(2 + _Z3 + 2 _Z4)

 _Z7 > -(1/4)*(2 _Z1 + 3 _Z2 + 2 _Z3 + _Z4 + 3 _Z5)

Since one now has lower bounds on _Z1.._Z8, and upper bounds are trivial, eight for -loops to test the inequality and all solutions can be generated. However is this noticeably "better" than running nine for loops for the variables in the original equation???

The only part of my original post which I *really*wish I had stated differently is the sentence

Would I guarantee that these are *all* the posint solutions - well no probably not

which I would now rewrite as

Would I guarantee that these are *all* the posint solutions - almost certainly not

In addition to the well known problems of local/global optimisation, I suspect that the DirectSearch package has an upper limit on the number of solutions it returns (although I cannot find this documented anywhere)

Theoretically, isolve() should provide the  answer which you want  - and it does if you remove the constraint that all variables are positive integers.

However even the constraint-free solution is given in terms of eight parameters, which is accurate, but probably does not represent significant progress on the original problem. I personally failed to make much fiurther progress with the isolve() command: not saying it can't be done, just observing that I got nowhere.

My next attempt was to use the DirectSearch package, which is not technically part of Maple, but is a package which you can download (from maple's website and install). A little experimentation with this and I came up with

   restart;
#
# rewrite the equation in a way which allows maxima and minima to be found (DirectSearch requirement)
#
   eqn1:=30*a+75*b+110*c+85*d+255*e+160*f+15*g+12*h+120*i-800000:
#
# Perform a global search which *should* return all combinations of posint variables which
# produce either maxima or minima
#
   ans:=DirectSearch[GlobalSearch](abs(eqn1), assume=posint):
#
# Just because I'm paranoid about max/min/opt routines, perform a quick dirty check
# on whether the returned solutions are valid (and are minima), return only those which are
#
   validAns:=[seq(`if`(eval(eqn1, ans[j][2])=0,ans[j][2],NULL), j=1..ArrayTools[Size](ans)[1])];
   numelems(validAns);

The result of this code was to produce 44 combinations {a..i} of positive integers which satisfy the original equation.

Would I guarantee that these are *all* the posint solutions - well no probably not. This is not a reflection on the author(s) of the DirectSearch package, just a reflection of my experience being often caught out by local/global optimisation problems.

It would be possible to check for further solutionsby feeding the DirectSearch ones to the parameterised isolve() solution mentioned above, to see if any more crop up

Checked this in Maple 18 and Maple2015, and it produces 4.0000000 in both.

Suggest

1. Put a restart command immediately prior to the g:=Groe.....  and check whether the error still occurs
2. Based on my tests the error will not occur - suggests that you have 'a' and/or 'b' defined to something a bit weird earlier in your worksheet

@jan123 

'sum' is for symbolic sums and it will do things like try to compute closed form solutions, (if it can) rather than just listing all the terms in the sum

'add' doesn't do any thing clever - it justs adds all the terms you give it

Respectfully suggest RTFM

@Thomas Richard 

As has been suggested

coeff(add(x^i,i=0..999)^500,x,2);

works (and keeps on working for more terms, higher powers). Didn't really attempt to find the limit doing it this way

Given Preben's comment I thought it might be possible to use the continuation method on the exp terms (see ?dsolve,numeric_bvp,advanced (bvp))

Having tried a couple of obvious "continuations", (including adjusting the continuation variable "manually") the *best* I managed to do was to convert the "Newton iteration is not converging" to "Error, (in dsolve/numeric/bvp) singularity encountered"

It may be possible to come with a continuation method which works, although I'm beginning to doubt it. I suggest you take a long hard look at the terms exp(beta/theta(eta)), particularly given the boundary condition theta(N) = 0

I still don't understand what you are getting trying to achieve: for example you say

"F might change depending on user input"

well F is local to your procedure p2() and values willl be assigned to it, depending on the y-values passed to p2()

Also you say that the "number of entries in p1 might change" What does this mean? That, as in your example p1() will always be passed a two element list but one might want to produce four(?) combinations of these two values- say

y[1]+y[2], y[1]-y[2], y[1]/2, y[2]/3 - easy to extend the computations in p1() to cover this case. Or do you mean that p1() will be passed an unspecified number of parameters and will contain and unspecified number of computations?? Very difficult to write a procedure which contains an unspecified number of calculations.

The best I can do is ask you to consider the following (about as close to your example as I can get) - copy/paste/run it and then tell me under what conditions it doen't do what you want

restart;
p1:= proc( y )
           return [y[1], y[2]-y[1]^2];
     end:
p2:= proc( y )
           local err, F;
           err:=5;
           if   y[1]<1
           then
              #
              # Call p1, and assign its output
              # to the local variable F
              #
                F:= p1(y);
           end;
         #
         # Print the local value of F just to
         # show that the correct assignment
         # has been made
         #
           printf("from print statement: F=%a\n",F);
         #
         # return the value of F which has been
         # computed - might not want to do this!
         # - might just want to use F further in
         # this procedure
         #
           return F;
     end proc:
#
xx:=p2([0.5, 7]);

@ben maths 

Because I have no idea what you are trying to achieve, I have no way of knowing why you want the the procedure root_of_cheb() as a sub-procedure or what argument you want to call it with!!!!!

If you have absolutely no idea when to call root_of_cheb() or with what argument - how I am expected to know???

@ben maths 

Following the guidelines I gave, you could write

restart:
root_of_cheb:=proc(n)
   local xk,b,k:
   xk:=(k,n)->cos((2*k-1)*Pi/(2*n));
   sort([seq(evalf(xk(k,n)),k=1..n)]):
end:

EvalInt:=proc(f,n)
                     local xk, root_of_cheb:
                     root_of_cheb:=proc(n)
                                           local xk,b,k:
                                                  xk:=(k,n)->cos((2*k-1)*Pi/(2*n));
                                                  sort([seq(evalf(xk(k,n)),k=1..n)]):
                                           end proc:
                    #
                    # Subsequent code code does not call root_of_cheb
                    # but it could, you just have to decide what argument
                    # to call it with!!!!!
                    #
                    xk:=(k,n)->cos((2*k-1)*Pi/(2*n));
                    evalf((Pi/n)*add(f(xk(i,n)),i=1..n)):
                    end:

Considered my previous reponse and decided that it was incomplete.

There are generally two scenarios for the operation you describe

Firstly a sub-procedure will be used multiple times within a given procedure, and will not be used anywhere else. In which case the sub-procedure can be defined as a local within the calling procedure: consider the example

restart;
dosqrt4x:=proc(x)
               local y, z, multBy2;
               multBy2:=proc(x)
                             return 2*x;
                        end proc;
               y:=  multBy2(x);
               z:=  multBy2(y);
               return evalf(sqrt(z));
          end proc:
dosqrt4x(5); #will compute sqrt(20)

Notice that multBy2() is defined as local within the procedure dosqrt4x() and is called twice. Since it is local to dosqrt4x(), it will not be avialbale to any other procedures

Secondly a sub-procedure is used in different procedures, in which case is necessary to consider it as a global: consider the example

restart;
multBy2:=proc(x)
              return 2*x;
         end proc:
dosqrt2x:=proc(x)
               local y;
               global multBy2;
               y:= multBy2(x);
               return evalf(sqrt(y));
          end proc:
dosqrt4x:=proc(x)
               local y;
               global multBy2;
               y:=2*multBy2(x);
               return evalf(sqrt(y));
          end proc:
dosqrt2x(5);#will compute sqrt(10)
dosqrt4x(5);#will compute sqrt(20)

Other combinations are available eg where a master procedure contains multiple sub-procedures and there is a utility procedure which is  required across the subprocedures, where the template would look something like

master:=proc()
             local utility, sub1, sub2;
             utility:=proc(...)
                        .......
                        end proc;
             sub1:=proc()
                         .......
                        #can use a call to utility()
                       end proc;
             sub1:=proc()
                         .......
                        #can use a call to utility()
                       end proc;
             #####
             code invlving calls to sub1(), sub(2) [or even utility)
             end proc;

Consider the following trivial example

multBy2:=proc(x)
              return 2*x;
         end proc;

dosqrt2x:=proc(x)
               local y;
               y:= multBy2(x);
               return evalf(sqrt(y));
          end proc;

dosqrt2x(10);

You need to set up a procedure which produces the correct answer irrespective of the vector which is supplied

So

z:=[a₀, a₁, a_2, a₁^*, a₂^*];
XYZ(z);

and

z:=[x,y,t,s,u]
XYZ(z)

are essentially equivalent - just a substitution of variable names

The attached is an attempt which *seems* to work. I have not attempted to simplify the code which produces the final matrix - yours seems extraordinarily lengthy/complicated for what it achieves!!!

I have added a simplified procedure at the end which does the same thing - I accept that this might/will fail if the supplied vector length is anything other than 5. It is intended to make you question whether all the code you have supplied is actually necessary.

atild2.mw

If I delete the option numeric from your pdsolve command then Maple produces an analytic solution. This would seem to be progress, but

1. The analytic solution contains four arbitrary constants, and since you only have three boundary conditions, substituting the boundary conditions into the analytic solution will still leave one condition short of a solution - four boundary conditions would seem to be necessary!
2. Even given four boundary conditions the analytic solution contains the LambertW function. I have no idea why this function crops up so often in Maple, but it is a complicated, multiple-valued function (check Wikipedia). I suspect, but cannot prove, that it occurs because of the existence of the quadratic term in diff(u(x,t),t) in your equation. Even given four boundary conditions, it may be that the multiple-valued nature of the LambertW function just confuses the hell out of the Maple numeric solvers, which I suspect are set up to deal with single-valued solutions. If so it suggest that it might be necessary to reconstruct the problem to avoid the possibility of multiple solutions. It is possible that ou *know* constraints on u(x,t) which might facilitate this, but since I do not, I cannot make much progress

@sami131 

If you had asked for a solution in Matlab, I would have been happy to provide it.

In fact you could have requested a solution in any of C/lisp/java/maple/matlab/perl/python (and a few others), and I would have considered it.

Trivial problems such as this can be solved in any programming language. Unfortunately sufficient  information is required for any of them - and any language would need the answers to my two original questions.

It is unfortunate (for you) that you cannot appreciate this simple fact