Add these options to the dsolve command:

method= bvp[midrich], maxmesh= 2^14

Also, here's an easier and much-more-accurate way to get your integral: Add the following ODE and BC to your system that you pass to dsolve:

diff(J(y), y)= E*diff(a(y), y),  J(-h/2) = 0

Then the integral is

1/2/(1+nu)*eval(J(y), aa(h/2));

and you can do away with the for loop and the PolynomialInterpolation.

If it's possible to solve at all, then I think that the following will do it:

P:= v-> cos(v)^2 + sin(v)^2 - 1:
V:= [Psi, Theta[1], Theta[2]](t):
U:= [(cos,sin)~(V)[], [x[1], x[2], z[1], z[2]](t)[]]:
thaw([eliminate(subs(U=~ freeze~(U), [[eqsys[], P~(V)[]], U])[])]);

This treats the sin and cos of Psi, Theta[1], and Theta[2] as separate variables and adds three equations based on the Pythagorean identity. So the system becomes a quadratic system of 10 equations in 10 variables. I let this run for over an hour and consume 5 Gig of memory, but I didn't have the patience to let it finish. I suggest replacing several of the parameters with small integers and trying again.

I used your ListParameters to set all parameters to 1, and I got no solution. So I think that there's a good chance that there's no solution.

There is a difference in the way Matrices evaluate between your Maple 15 and the current Maple 2016. I don't know the precise details, but you need to add something to VV's command in order for it to work in your Maple. Change

C[p0]:= eval(C[iso], const);

to

C[p0]:= eval(rtable_eval(C[iso]), const);

I have verified that

• the command that VV gave works in Maple 2016. That is, it produces a fully evaluated numeric Matrix;
• VV's command by itself doesn't work in Maple 16;
• and the modified command that I gave works in Maple 16.

It's very difficult for me to describe to you what exactly rtable_eval does. It's a bit mind boggling.

You wrote:

The resulotion of plot is not so good, and when I export this plot as jpeg file, quality of resolution deteriotate even more, I mean pixel tear out on zooming the image. How I fix this issue which I guess might be due to wrong way of exporting plot ?

JPEG is not a good format for most plots. JPEG is intended for photos (the P in JPEG stands for Photographic), which have continuous color changes. Plots tend to have crisp, abrupt, discontinuous color changes, especially in the transition from foreground to background. The "pixels tearing out" is JPEG's attempt to convert those abrupt changes into continuous ones.

On the other hand, the G in both PNG and GIF stands for Graphics. These are better than JPEG for most plots.

You had many syntax errors---too many for me to list right now---and your procedure CreaF would've been very inefficient if it had worked. Here's a better way:

f:= (x,y,i,j)-> ((x+j-1)/2^i, (y+2^i-j-1)/2^i):
CreaF:= (i::posint)-> unapply(Vector(2^i-1, j-> f(x,y,i,j)), (x,y)):

And call it like this:

CreaF(2)(.35,.465)[1];

     0.875000000000000e-1, .616250000000000

Note that I changed the call order from your CreaF(2)[1](.35,.465) because I thought that that would generally be more efficient. In other words, CreaF now returns a procedure which returns a Vector as opposed to CreaF returning a Vector of procedures. If you need it to be the latter way, let me know.

Your implementation of Kitonum's Answer is much longer than it needs to be, and takes much more time. Here's my one-liner procedure for it:

DeleteZeroColumns:= (M::Matrix)-> M[.., remove(k-> M[1,k]=0, [$1..op([1,2],M)])]:

This runs in 3 to 4 milliseconds on an 8x3000 Matrix. The code op([1,2],M) returns the number of columns of M.

If the first argument to display is an array of plots, then you'll get a table of plots. If the first argument is a list of plots, then you'll get them all plotted in the same coordinate plane.

In Maple 2016: Suppose that A is the 5x2 Matrix. Then do

[seq(convert(v[], list), v= Iterator:-CartesianProduct(convert(A, listlist)[]))];

To see the full solution for 2nd-4th degree polynomials with variable coefficients, use option explicit:

solve(myeq, y, explicit);

Use PDEtools:-dchange to convert to polar coordinates, where the boundary condition is simpler.

PDE:= diff(u(x,y), x$2) + diff(u(x,y), y$2) = 1:
PDEp:= PDEtools:-dchange({x= r*cos(theta), y= r*sin(theta)}, PDE, simplify);

     PDEp := ((diff(u(r, theta), r, r))*r^2+(diff(u(r, theta), r))*r+diff(u(r, theta), theta, theta))/r^2

pdsolve({PDEp, u(1, theta) = 0});

     u(r, theta) = (1/4)*r^2+_F1(I*ln(r)+theta)+_F2(-I*ln(r)+theta)

Change y^`*` to `y*`.

(Some of this will repeat what Joe said, but I have the improved code also.)

Your hunch that hash tables would be better is spot on. Although your algorithm is impressive, your Maple implementation of it suffers from the most common of newbie efficiency mistakes: building a set or list by iteratively appending to an existing set or list. (I even used to do it myself, many, many years ago.) This can be corrected by building a table instead and then using indices or entries to build the set or list.

It only required a trivial change to make this change for your big set, orbits. This simple change reduced the run time to under 15 seconds. By also correcting the same mistake for your small set orbit and converting your small list interl to one built with a single seq, I shaved off another three seconds. By converting some of the other for loops to seqs, I shaved off another second. Moving the computation of d2 outside the inner loop saved another three seconds (an idea that I got from Joe). Remembering the previously computed values of d3 saved another 6 seconds (an idea that I got from Joe). The final code's memory usage is 22.18 M, and its time is under 3 seconds. Here is the code:

HEX:= proc()
local orbit, orbits:= table(), idx2, idx3, d2, d3, pos,
    interl, r, entry, refl, q,
    Base3:= proc(n) option remember; convert(n, base, 3) end proc
;
    for idx2 from 2^6 to 2*2^6-1 do
        d2:= convert(idx2, base, 2);
        for idx3 from 3^6 to 2*3^6-1 do
            d3:= Base3(idx3);
            interl:= [seq([d2[pos], d3[pos]][], pos= 1..6)];
            orbit:= table();
            orbit[seq(interl[[seq(q, q= 1+r..12), seq(q, q= 1..r)]], r= 0..11, 2)]:= ();
            for refl from 0 to 11 by 4 do
                entry:= interl[[seq(q, q= 1+refl..12), seq(q, q= 1..refl)]];
                orbit[entry[[1, seq(12-q, q= 0..10)]], entry[[3, 2, 1, seq(12-q, q= 0..8)]]]:= ()
            od;
            orbits[{indices(orbit, nolist)}]:= ()
        od
    od;
    nops([indices(orbits, nolist)])
end proc:

This code isn't intended to be run as 2D input, which fails by prematurely evaluating ().

You already have the correct implicit solution. It's what you called implicitsoln1. There's no way to get an explicit solution, i.e., a solution without RootOf and with y isolated on one side of the equation, which is probably why your instructor asked you for the implicit solution.

In the expression D(u)(x,0), you need to specify the independent variable with respect to which the derivative is being taken. So, you can change it to D[1](u)(x,0) for the derivative with respect to x or D[2](u)(x,0) for the derivative with respect to t.

You can start by reading ?HelpTools and its subpages. You should create your help file as a worksheet rather than as a .txt file.