@Mac Dude 

To avoid the interpolation one may take the inverse transform (provided it is bijective).
So, let A be a subset of the rectangle R(A), and f:A-->B  the transform (bijective).
B subset of a rectangle R(B).
Simply take the color(b) := color(g(b)), b in B, where g is the inverse mapping.

@Markiyan. Unortunately plottools:-transform can be applied only to the plot structure, not to the background.
One may try plot3d(...,image=filename).

@Mac Dude 

Conformal mappings are easy to construct:

(x,y) |--> (u,v) = (Re f(x+iy), Im f(x+iy)),

f being holomorphic such as f(z)=z^2, f(z)=1/z, f(z)=sin(z) etc.

[ The example is essentially (x,y) |--> (x, y^3) which is of course not conformal;
a conformal example would be e.g.  (x, y) |-->  (x^2-y^2, 2*x*y) ]
The origin does not matter, a translation does not affect the conformality.

Of course the image will have to be clipped (into a rectangle); if f(z_0) = oo, it wil be of course clipped.
The example can be easily adapted.

Thanks, I have inserted.

@Preben Alsholm 

@rlopez 

I think that the (global) bivariate limits for rational functions were introduced in Maple 2015 only.
I remember this because I was really impressed by the new algorithm.

limit( f(x,y), {x=0,y=0}); # (1)

is the limit (multi) at (0,0).
If it exists, it equals

limit( f(r*cos(t), r*sin(t)), r=0); # (2)

But it is possible that (2) exists and (1) does not (because (2) is not necessarily uniform wrt t).
For example, consider f(x,y) = x^2/y  or  f(x,y) = x^2*y/(x^4+y^2).

@Preben Alsholm 

I would vote for it, but I am not affiliated to Maplesoft :-)

Actually Maple uses residues in computing integrals, see e.g.

int(ln(x)/(x^6-2*x+4), x=0..infinity);

@Marko Riedel 

In this case it seems that using residues is the only solution:

where L = log, R = rational function without zeros in [0,infinity)
(this is what you did, J(n) being so expressed via J(n-1), J(n-2),... ).

@Marko Riedel 

tedious to code?

pw:=(n,a) -> add(NumberTheory:-Moebius(n/d)*a^d, d in NumberTheory:-Divisors(n)):

Test:
StringTools:-Generate(4, "01"):
select( StringTools:-IsPrimitive, %);
pw(4,2)=nops(%);

@MortenZdk 

The fact that your command did not work is more a Maple's fault than yours. The Explore command is rather new and must still be polished. In Maple 2015 not even the f:=a*x ...  version did not work, it had to be entered as

Explore(plot(a*x, ...

So, keep studying and try the examples in the Help (you can turn them into worksheets so that you can experiment). Also, User's Manual is a good start.

Of course it is possible. But:

1. This is equivalent to det(J)=0
2. In solving such systems, some solutions will be lost.
    See e.g.  solve(x*y+z*t);

Since you did not provide an example, I am not able to help more.

@cuongtd 

@Preben Alsholm OK, it's collocation then :-)
But usually in collocation methods the number of points equals the number of coefficients and a system is solved.
In my case I took nx >> n, that's why I think that choosing "special" points should not matter.
(I was not concerned about efficiency).

@Preben Alsholm 

With fsolve, the situation is clear.
The problem is with the minization approach.
With your example, this approach should find the exact solution (being in the basis). Here it's only NLP's fault of not being able to minimize properly.

I have tested also another example:

ide:=x^2*(diff(y(x), x, x))+(71/12)*x*(diff(y(x), x))-(7/4)*y(x)+int(y(t), t = 0 .. 1);
exactsol:=x^(1/3)+1;

which is not in the basis x^(k/2).

Added. By "no collocation" I meant not special chosen points.

@Preben Alsholm 

I tested you example withot collocation points and I got |y - yexact| < 2*10^(-6)

(This is normal because the exact solution is - unfortunately - in the basis!).

I wonder about the benefit of using collocation points instead of equidistant ones in this approach
(because the minimization is used).

@Preben Alsholm 

I just wanted to do something similar.
An analytic example would be nice.

Added. If someone is interested, here is y2 (without Chebychev polynomials and collocation)

Digits:=15:
ide:=x^2*(diff(y(x), x, x))+50*x*(diff(y(x), x))-35*y(x) - (1-exp(x+1))/(x+1)-(x^2+50*x-35)*exp(x)-int(exp(x*t)*y(t),t=0..1):
n:=10: nx:=16:
p:=x->add( a[k]*x^(k/2),k=0..n): 
F:=unapply( eval(ide, [y=p]),x):
d:=evalf(add(F(k/nx)^2,k=1..nx)):
s:=Optimization:-NLPSolve(d, {p(0)=1, p(1)=1.5}, iterationlimit=20000):
yy:=eval(p(x),s[2]);   # approx solution
yy=
   1.+.128100672246220*sqrt(x)-4.96593663847107*x+21.9900671494260*x^(3/2)-
   61.2588096618634*x^2+119.676248631378*x^(5/2)-158.798014306812*x^3+
   142.246959241914*x^(7/2)-81.9327391785196*x^4+27.4658220286035*x^(9/2)-4.05169793790087*x^5
The plot is the same as Preben's.