Here's an example of a procedure that uses the parse statement to create another procedure with for loops nested to arbitrary depth. The goal of the latter procedure is to generate the Cartesian product of n lists, n being the argument of the first procedure. This is metaprogramming: writing a program that writes another program. It is a profoundly powerful technique which is much easier in Maple than in most other languages.

CartProdGen:= proc(n::posint)
local d; #depth
     parse(
          cat(
               "proc(", seq(cat("L",d,","), d= 1..n), "$)",
                    "local n:=0,",seq(cat("i",d,","), d= 1..n),
                    "R:= table()",
               ";",
                    seq(cat("for i", d, " in L", d, " do "), d= 1..n),               
                         "n:= n+1;",
                         "R[n]:= [", seq(cat("i",d,","), d= 1..n-1), "i", n, "]",
                    "end do " $ n, ";",
                    "convert(R, list)",
               "end proc;"
          )
     )
end proc:

Example of use:

CartProd3:= CartProdGen(3); #Generate the procedure

CartProd3([a,b], [c,d], [e,f]); #Use it

The two steps can be put together into a single procedure so that the generated procedure is immediately used:

CartProd:= (L::list(list))-> CartProdGen(nops(L))(L[]):

CartProd([[a,b], [1,2], [c,d]]);

There are many completely different very short procedures (several as short as two lines) for the Cartesian product of an arbitrary number of lists that don't use any metaprogramming (don't use the parse statement). They can be found by doing a MaplePrimes search on "Cartesian product of lists".

This is a very easy operation in Maple. Example:

R:= rand(0..9):
L:= ['R'() $ 100]:
nops(L);

     100

Unique:= {L[]};

CountUnique:= (L::list)-> nops({L[]}):
CountUnique(L);

     10

The primary way that a procedure returns a value is through a return statement, not through a variable declared evaln. The latter way is an obscure technique used for procedures that need to return something in addition to their primary return value. Indeed, it's hardly ever used now that procedures can return an expression sequence.

Max:= proc(L::list)
local i, Max:= -infinity;
     for i to nops(L) do
          if L[i] > Max then
               Max:= L[i]
          end if
     end do;
     return Max
end proc:

If the return statement is the last statement, then the word return can be omitted. So the above could simply have Max on the second-to-last line.

There's no syntactic reason why I used Max as both a local variable and a procedure name; it's just my style.

I avoided using max as the procedure name because it would conflict with Maple's built-in max. Overwriting the name of a built-in procedure is allowed, but that's another topic.

In Maple, Pi is spelled with a capital P. Your cotx will need to be changed to cot(x). You will need to supply numeric values for D1, D2, a, and q before you can get a solution. The highest derivative of P in your ODEs is the first. Thus, you can only have three initial conditions, and you can't use D(P) as one of them.

Your particular solution, (i.e., your solution to the inhomogenoeus part) is correct. To get the general solution, you'll need inital conditions. In the code below I've used initial conditions that y, y', y'', and y''' are all 0 at t=0.

u:= exp(3*t)+3*exp(t):
ode:= diff(y(t),t$4) - 16*y(t) = diff(u,t) + u:
ICs:= seq((D@@k)(y)(0) = 0, k= 0..3):

dsolve({ode, ICs});

ee:= sin(3*arcsin(x)):
expand(%);

Use command Bits:-Join. This will require you to enter the number as a list of 0s and 1s. If you'd rather not enter the commas that will be required for that, here is a very small procedure that will convert the input from a string of 0s and 1s.

BitStringInput:= Bits:-Join@parse~@StringTools:-Explode@StringTools:-Reverse:
x:= BitStringInput("10010");

 

You can perform bitwise logical operations on your bit strings (which, for efficiency, are represented internally as integers) using other commands in the Bits package. When you want to see your output as a bit string, do

Out:= Bits:-String(x, msbfirst);

restart:
f:= (x^2+2*x+3) - (k+5*x-7*x^2):
solve(int(f, x= `..`(solve(f, x))) = 36, k);

Here is a way completely different from Kitonum's. In the procedure below, I take the two quadratic solutions returned by solve and determine which is the positive branch numeric evaluation.

PositiveRoot:= proc(
     eqn::{algebraic,`=`},
     var::name,
     {evals::set(name= complexcons):= indets(eqn, And(name, Not(constant))) =~ 0}
)
local
     Sol:= [solve(eqn, var)],
     E:= evalf(eval(Sols, evals))
;
     if nops(E) = 0 then [][]
     elif nops(E) = 1 then Sol[]
     elif nops(E) > 2 then
          error "Procedure only intended for quadratic equations."
     elif Im(E[1]) > Im(E[2]) then Sol[1]
     elif Im(E[1]) = 0 then
          `if`(E[1] > E[2], Sol[1], Sol[2])
     else
          Sol[2]
     end if
end proc:

Example of use:

PositiveRoot(eq1, R);

I echo Markiyan's concerns: The extremely large coefficients, especially in eq3 and eq4, overwhelm the significance of the relatively small final results. Consider using exact arithmetic.

The proximate cause of your error message is that you enclosed the boundary conditions in square brackets. Thus, you have two sets of brackets in the call to pdsolve. 

Any finite, simple, closed region of a computer plot is a polygon, being composed of a finite number of line segments, however small. Thus it can be filled with plots:-polygonplot. You may need to extract the vertices from the original plot structure to do this. If you post your code, I can show you how. 

Indeed A^%T is not even necessarily square, let alone invertible. Use command LinearAlgebra:-LeastSquares(A, b, optimize). You can also multiply b by the pseudoinverse of A, which you can obtain as LinearAlgebra:-MatrixInverse(A, method= pseudo).

Your syntax is perfect, and it gives me the answer -3/2. Try using the restart command.

Taking the constants to be 1 and plotting shows that there is only one positive solution. The solve gives me one answer plus the warning. So, just ignore the warning and use the answer, which is h*c/T/k/(LambertW(-5*exp(-5) )+5).

I'm not familiar with the coeffs option to dsolve. Where is it documented? If I try it, I just get an error message that leads me to believe that no such option exists. But just trying without any options, I get an answer in terms of integrals of Airy functions.