Avoiding automatic simplifications is a very tricky thing, and I don't think that there's any perfect solution. Your case is especially tricky because you also want to preserve superfluous parentheses. I think that the best solution is offered by the InertForm package.

macro(IF= InertForm):

Enter the expression as below. The purpose of the ``s is to preserve the parentheses.

ex:=
     ``(IF:-Parse("F^5*alpha[5]+F^4*alpha[4]+F^3*alpha[3]+F^2*alpha[2]+F*alpha[1]")) +
     ``(IF:-Parse(
            "F^4*G*gamma[1]+F*G^4*gamma[3]+F^2*G*gamma[2]+"
            "F*G^2*gamma[4]+F*G*gamma[5]"
     )) +
     ``(IF:-Parse("G^5*beta[5]+G^4*beta[4]+G^3*beta[3]+G^2*beta[2]+G*beta[1]")) +  
     delta[0]
:

To view the expression, use

IF:-Typeset(ex, inert= false);

To use the expression mathematically, use

expand(value(ex));

It's not the complexity that matters. The problem is that (as far as I can tell) geometry:-draw can only plot within the window [-10..10, -10..10]. Your first line doesn't intersect that window.

As a workaround, use implicitplot:

with(geometry):
point(P1, [47+(38+22/60)/60, -(122+(43+4/60)/60)]):
point(P2, coordinates(P1) +~ [cos,sin](30*Pi/180)):
#Specify axes names as 3rd argument to `line` to avoid super-annoying dialogs!
line(L1, [P1,P2], [x,y]):
plots:-implicitplot(Equation(L1), x= -100..100, y= -100..100);

Update: The second sentence above is incorrect: The default view is [-10..10, -10..10], but that can be changed or nullified.

Change the line D:= A to D:= copy(A). The line D:= A doesn't mean what you think it does. It means that D will be A and will stay A until it is reassigned. Thus, any changes to D also change A.

There's a problem with how you entered T: You can't use square brackets for algebraic groupings. Enter T like this:

T:= (r,theta)-> 0.5*alpha+Sum(gamma*n*(r^n+R^(2*n)+r^(-n))*sin(n*theta), n= 1..infinity):

Use PDEtools:-dchange to make the change of variables:

PDEtools:-dchange({theta= arctan(y,x), r= sqrt(x^2+y^2)}, Diff(T(r,theta), y), [x,y]):
value(%);

Notice that Diff is with a capital D. This is necessary to delay the evaluation of the derivative until after the change of variables is applied. The value then performs that evaluation, which takes about three minutes, unfortunately.

What you call a perfect matching, I call a bijection. Here's a little procedure that finds all bijections between lists A and B:

Bijections:= proc(A::list, B::depends(And(list, satisfies(B-> nops(B)=nops(A)))))
     map[3](zip, `[]`, A, combinat:-permute(B))
end proc:

Bijections([1,2,3], [4,5,6]);

[[[1, 4], [2, 5], [3, 6]], [[1, 4], [2, 6], [3, 5]], [[1, 5], [2, 4], [3, 6]], [[1, 5], [2, 6], [3, 4]], [[1, 6], [2, 4], [3, 5]], [[1, 6], [2, 5], [3, 4]]]

If you have Maple 2016, then the following may be more efficient:

Bijections:= proc(A::list, B::depends(And(list, satisfies(B-> nops(B)=nops(A)))))
local b;
     seq(zip(`[]`, A, convert(b[],list)), b= Iterator:-Permute(B))
end proc:

Remove the line

xi:= k[1]*x-k[1]*c[1]*t ;  eta:= k[2]*x-k[2]*c[2]*t ;

Replace it with

PDEtools:-dchange({xi= k[1]*x-k[1]*c[1]*t, eta= k[2]*x-k[2]*c[2]*t}, Diff(u,t), [x,t]):
value(%);

In Maple, the end of a line is treated like a space, so there's no need for any line-continuation symbol. For human readability, the continued parts of the lines should be indented, but the code parser doesn't care. In Maple, your example command could be

A:= [
     3, 4, 5, 6, 6, 45, 37,
     5, 4, 67, 39, -967
];

In the extremely rare situation that you need to end a line at a place where a space wouldn't be allowed, end it with a backslash \. But I can't think of any situation where there's not a more-elegant way to do this.

I assume that you mean that you want to import the drawing of the graph rather than its mathematical description.

First check which formats the other program can import. Left-click on the drawing so that it's selected, but so that no subpart of it is selected (a selected subpart would be highlighted). Right-click to bring up the context menu. Select Export, and select the format from the submenu. If possible, avoid JPEG: it's not a good format for line drawings.

It's possible to automate this process so that it's done programmatically. If that's what you want, let me know.

< <Matrix(S*(dmax-1), S) | Matrix(S*(dmax-1), shape= identity)>,
    -`<|>`(CCnew||(0..dmax-1))
>;

Explanation:

<...|...> or `<|>`(...) is used to separate columns or to separate blocks horizontally; <...,...> or `<,>`(...) is used to separate rows or to separate blocks vertically.

A sequence of global variables with numeric suffixes can be created as V||(a..b).

For reasons that I don't understand, the minus sign can't be right next to CCnew.

The assumptions unfortunately don't help in this case. The equations and inequalities need to be handled separately (as far as I know). Use command isolve on the equations. The Z is to parametrize the solutions.

Sol:= isolve({eq1, eq2}, Z);

     Sol:= {A = 39-49*Z, J = 60+40*Z, T = 1+9*Z}

The parameter Z can be any integer. We need to get the range of Z that satisfies the inequalities.

SolZ:= isolve(rhs~(Sol) >=~ 0);

     {Z = 0}

eval(Sol,  SolZ);

     {A = 39, J = 60, T = 1}

spin:= (-1)^~LinearAlgebra:-RandomMatrix(N, generator= 0..1);
to 1000 do
     neighbors:= ArrayTools:-CircularShift(spin, 0, 1) + ...
     ...
end do;

Is this what you want?

u:= x(t)/sqrt(x(t)^2+y(t)^2-2*x(t)+1)*diff(y(t),t);
v:= subsindets(u, specfunc(anything, diff), freeze);
T:= indets(v, name);
w:= subsindets(v, function(identical(T[])), freeze);
var:= indets(w, name);
mtaylor(w, var, 4);
thaw(%);

I got it to work using Groebner:-Basis instead of Primfield. I assume that you still want the fourth root of the starting matrix; that's what's computed in the worksheet below. Other roots could be computed just as well.

I couldn't have done this without the crucial idea supplied by Christian Wolinski, so I promoted his Reply with that idea to an Answer so that I could give it a vote up.

Sorry, for some reason this worksheet won't display inline, but here's the link:

Matrix_powers_finite_field_2.mw

It seems not possible to get a solution for some negative values of T. The following code will produce a matrix of the solutions that it is possible to get.

restart:
Digits:= 15:

eq1 := PI_T = alpha*M*(kappa/(1+kappa-sigma)-1):
eq2 := l = alpha*((sigma-1)*kappa/(1+kappa-sigma)+1):
eq3 := M = phi_c^(-kappa)*F:
eq4 := e*F = PI_T*v+T:
eq5 := M*l+L_A = L_s:
eq6 := F*e+A = A_s:
eq7 := A = (1-beta)*(L_s+(1-v)*PI_T-T):
eq8 := A_s = L_A:
eq9 := p = sigma/((sigma-1)*phi):
eq10 := P_Y = M^(lambda/(1-sigma))*p:
eq11 := L_s = N*(theta*P_Y^beta)^(-1/delta):
eq12 := phi = (kappa/(1+kappa-sigma))^(1/(sigma-1))*phi_c:
eq13 := U = delta*theta*L_s/((1+delta)*P_Y^beta)+((1-v)*PI_T-T)/P_Y^beta:

Params:= {N = 1, v = 1, beta = .75, lambda = 1, sigma = 5, kappa = 4.8, delta = 2, theta = 2.591350635, alpha = 0.3998699153e-1, e = 0.8333333332e-1}:

Init_Values:= {A_s = .2500000000, l = .9996747882, PI_T = 0.8333333332e-1, L_A = .2500000000, phi = 1.878101496, M = .4168022157, A = .1666666667, p = .6655657336, P_Y = .8283396488, L_s = 2/3, phi_c = 1.2, F = 1, U = 1.326439132 }:

V:= <T, A_s, l, PI_T, L_A, phi, M, A, p, P_Y, L_s, phi_c, F, U>^+:

SOLs:= table():
for _T from -0.3 by 0.01 to 0.3 do     
     s:= fsolve(eval({eq||(1..13)}, Params union {T = _T}), Init_Values);
     if op(0, eval(s, 1)) = fsolve then next end if;
     SOLs[_T]:= eval(V, s union {T= _T})
end do:

sol_matrix:= Matrix(<entries(SOLs, 'nolist', 'indexorder')>, datatype= float[8]);

simplify(P(epsilon), {epsilon^2= 0});

This will change epsilon^2 and all higher powers of epsilon to 0.