The expression RootOf(expr) stands for the root(s) of expr. The _Z appearing is the unknown in expr.

Try entering the following and see what you get.
RootOf(x^2-1,x);

The 3-argument version of signum is explained in the help page for signum in the 3rd paragraph in the Description section.

To get numerical values use evalf.

To avoid the appearance of RootOf in this and similar cases you can add the argument 'explicit' to solve:

restart;
X:=l1*cos(theta1)+l2*cos(theta1+theta2)=x;
Y:=l1*sin(theta1)+l2*sin(theta1+theta2)=y;
res:=solve({X,Y},{theta1,theta2},explicit);
params:={l1=5,l2=5,x=5,y=3};
eval(res,params);
evalf(%);
##Now without 'explicit'
res2:=solve({X,Y},{theta1,theta2});
eval(res2,params);
evalf(%);

convert(BesselJ(0,x),Sum);
subs(_k1=n,convert(%,factorial));

Well, there is at least one more. What is your purpose in doing this?

Try this.

restart;
eq := diff(f(eta), eta, eta)+(1/(4*eta^2)+(3*n+1)/(n+1)*(NU/(4*a*eta)-NU*(eta/a)^((n+1)/(2*n))/(4*a*eta)))*f(eta) = 0;
dsolve(eq);
#Returns a DEsol structure, which can be used for some purposes, see the help
#Why did you get an answer for n=1/5?
eval(eq,n=1/5);
dsolve(eval(eq,n=1/5));
#Maybe because the exponent is an integer. So we try that:
solve(p=(n+1)/n/2,n);

##We must exclude p=0 as then n+1=0

for p in [seq(-5..-1),seq(1..5)] do res[p]:=dsolve(eval(eq,n=1/(2*p-1))) end do;

Success for p=-2, -1, 1, and 3. The result for p=3 (i.e n=1/5) is not quite explicit since it still has an unevaluated integral.

Replace sign with signum.

Then just do as you did:

y:=0: eq:=0:
for j from 1 to ny do
  y:=(j-1)*delt_y+delt_y/2;
  eq:=evalf(eq+b(y)*signum(y-y0)*abs((y-y0))^(1/n)*delt_y)
end do:
sol1:=fsolve(eq=0,y0); #Or the syntax you used

For b the simpler definition

b:=y->piecewise(y<=H1,b1,y<H1+H2,b2,y<=H,b1);

accomplishes the same as yours assuming that H2>0 and H >=H1+H2.

Adding a view option helps:

plot(1/cos(t), t=0..Pi/3, coords=polar,view=[0..2,0..2]);

It seems that the view used by default is determined by the max and min of the x- and y- values:

plot(1/cos(t), t=0..Pi/3, coords=polar);
plottools:-getdata(%);
itv,M:=%[2..3][];
plot(M,view=itv);

(max,min)(M[..,1]);
(max,min)(M[..,2]);

I think that the problem is numerical and cannot be avoided.
Consider a ball bouncing just up and down. Assume that on hitting the ground the speed is reduced by a factor k < 1. Then the total time T for infinitely many bounces is finite. Thus it is not surprising that any numerical procedure is going to fail eventually and presumably before time T has passed.

restart;
ode:= diff(y(t), t, t) = -g;
ics := D(y)(t0) = y1,y(t0) = 0;
sol := dsolve({ode, ics},y(t));
solve(rhs(sol)=0,t);
#Time spent in the air is proportional to the starting velocity.
#Upon hitting the ground the starting velocity is multiplied by k (k<1).
#Thus the time it takes for N bounces is
Sum(2*y1*k^i/g,i=0..N-1);
#Notice that the time for infinitely many bounces is finite:
sum(2*y1*k^i/g,i=0..infinity) assuming abs(k)<1;
#This is going to fail at some point
sol:=dsolve({eval(ode,g=1),y(0)=0,D(y)(0)=1},numeric,events=[[y(t)=0,diff(y(t),t)=-diff(y(t),t)*0.8]]);
plots:-odeplot(sol,[t,y(t)],0..11);
#The finite time for the infinitely many bounces is in this case:
sum(2*1*0.8^i,i=0..infinity); #returns 10

########### In this simple case where T is easily calculated you could just add an event:

sol:=dsolve({eval(ode,g=1),y(0)=0,D(y)(0)=1},numeric,events=[[y(t)=0,diff(y(t),t)=-diff(y(t),t)*0.8],[t=10,halt]]);
plots:-odeplot(sol,[t,y(t)],0..11);

_z1 is just an integration variable, just as t is in this integral:

int(t^2, t=0..x);

The reason for the somewhat strange name is that dsolve has to find a global name itself, one which doesn't conflict with names you are using or would be likely to be using.

You can try value(f) to see if Maple can evaluate the integral: Notice that it is grey, which indicates that it represents Int, the inert version of int.

It is obviously much better to get real code instead of pictures, but you seem to be aware of that.

Using events with (essentially) your piecewise version although via a discrete variable b(r):

restart;
polytrope1:=4./3;
polytrope2:=2.;
K2:=10^12;
rho_crit:=4.*10^(-14);
K1:=K2*rho_crit^(polytrope2-polytrope1);
odes := diff(rho(r), r) = -(rho(r)*(1+K*rho(r)^(polytrope-1)/(polytrope-1))+K*rho(r)^polytrope)*(4*Pi*r^3*K*rho(r)^polytrope+m(r))/(K*polytrope*rho(r)^(polytrope-1)*r*(r-2*m(r))), diff(m(r), r) = 4*Pi*rho(r)*(1+K*rho(r)^(polytrope-1)/(polytrope-1))*r^2;
r0:=10.^(-10);
ics:=rho(r0)=10.^(-13),m(r0)=0;
frho1,fm1:=op(subs(K=K1,polytrope=polytrope1,rhs~([odes])));
frho2,fm2:=op(subs(K=K2,polytrope=polytrope2,rhs~([odes])));
frho:=b(r)*frho1+(1-b(r))*frho2;
fm:=b(r)*fm1+(1-b(r))*fm2;
sys:=diff(rho(r),r)=frho, diff(m(r),r)=fm;
ics2:=ics,b(r0)=1;
##with b(r0)=1 we start with K1 and polytrope1. That is changed when the event rho(r)=rho_crit happens to b(r)=0.

res:=dsolve({sys,ics2},numeric,discrete_variables=[b(r)::boolean],events=[[rho(r)=rho_crit,b(r)=0]]);
1e16*r0;
plots:-odeplot(res,[r,rho(r)],r0..1e6);
plots:-odeplot(res,[r,m(r)],r0..1e6);
plots:-odeplot(res,[r,b(r)],r0..1e6);

I tested the code in Maple 18 and Maple 15. It worked just fine in both.

I don't have access to Maple 11 at this machine, but do have access to Maple 12.
The code didn't work, but came up with an error message about being
"unable to handle PDE problem subject to boundary conditions {u(x,0)=phi(x)} "

I have no access to Maple 13 or 14.

I don't know why it fails to find a solution, but you can obviously do it yourself in steps, here aided by Maple:

restart;
PDE:=diff(u(x,y),x,y)=sin(x)*cos(y);
pdsolve({PDE,D[1](u)(x,Pi/2)=2*x,u(Pi,y)=2*sin(y)}); #No solution found
sol:=pdsolve(PDE); #The general solution
eq1:=eval(rhs(sol),x=Pi)=2*sin(y);
eq2:=eval(diff(rhs(sol),x),y=Pi/2)=2*x;
res1:=solve(eq1,{_F1(y)});
res2:=dsolve(eq2);
res3:=eval(res2,x=Pi);
res:=subs(res1,res2,res3,sol);
## Checking:
pdetest(res,PDE);
eval(res,x=Pi);
eval(diff(res,x),y=Pi/2);

Is this what you had in mind?

restart;
Xsubs:=X=x*exp(epsilon*alpha[1]);
Ysubs:=Y=y*exp(epsilon*alpha[2]);
psubs := P(X,Y) = p(x,y)*exp(epsilon*alpha[3]);

Eq:= diff(p(x,y),y$2)+diff(p(x,y),x)=0;

psubs2:=op(solve(psubs,{p(x,y)}));
Xsubs2:=op(solve(Xsubs,{x}));
Ysubs2:=op(solve(Ysubs,{y}));

PDEtools:-dchange({Xsubs2,Ysubs2,psubs2},Eq,[X,Y,P]);
combine(%);

####  Or are the parameters alpha and epsilon known functions of x and/or y ?

Use Optimization:-LSSolve on the residuals, i.e. the differences between the left and right hand sides of your g's.

restart;
Digits:=15:
Rx := Matrix (3,3, [1,0,0,0,cos(a),-sin(a),0,sin(a),cos(a)]);    
Ry := Matrix (3,3, [sin(b), 0, cos(b), 0, 1, 0, -sin(b), 0, cos(b)]);
Rz := Matrix (3,3, [cos(c), -sin(c), 0, sin(c), cos(c), 0, 0, 0, 1]);
RT := Rx.Rz.Ry;
#You didn't provide mat so I'm making it up:
mat:=eval(RT,{a=.123,b=1.5,c=-2.});
resids:=[seq(seq(RT[i,j]-mat[i,j],i=1..3),j=1..3)]; #The residuals
res:=Optimization:-LSSolve(resids);
#Check
eval(RT,res[2])-mat;

After experiments I'm guessing that the answer for an nxn- matrix is:

If n is even then all eigenvalues are real for abs(x)<=1. Outside that x-interval there are n-2 real eigenvalues only.

If n is odd ( = 2*k+1 ) then all eigenvalues are real for abs(x)<= sqrt((k+1)/k). Outside that x-interval only n-2 eigenvalues are real.

Can you prove that?

In experimenting try the code below for lower values of n to gain insight.

restart;
Digits:=30; #Necessary for fsolve for high values of n.
A:=n->Matrix(n,{(1,1)=-I*x,(n,n)=I*x, seq((i,i-1)=-1, i=2..n), seq((i-1,i)=-1, i=2..n)}):
p:=n->sort(LinearAlgebra:-CharacteristicPolynomial(A(n),lambda));
plots:-implicitplot(p(10),x=-5..5,lambda=-5..5,gridrefine=2,grid=[50,50]);
#The guess is primarily based on the existence of the closed curve in the plot made above.
#The plot is shown at the bottom.
#The guess is that the crucial dividing x-points are the values of x for which lambda=0.
q:=n->collect(p(n),x,factor);
LE:=[seq(eval(q(n),lambda=0),n=2..20,2)];
LO:=[seq(eval(factor(q(n))/lambda,lambda=0),n=3..21,2)];
solve~(LO,x);
q(50);
eval(%,lambda=0);
fsolve(eval(q(50),x=1.0001),lambda);
nops([%]);
fsolve(eval(q(50),x=0.9999),lambda);
nops([%]);
fsolve(eval(q(50),x=1),lambda);
nops([%]);
#49=2*24+1;
fsolve(eval(q(49),x=sqrt(25/24)+0.0001),lambda);
nops([%]);
fsolve(eval(q(49),x=sqrt(25/24)-0.0001),lambda);
nops([%]);
fsolve(eval(q(49),x=sqrt(25/24)),lambda);
nops([%]);

If the resulting ode is not required to be autonomous, i.e. the right hand side allows the appearance of explicit x-dependence as in your example, then the question is ridiculous or as Rouben Rostamian more politely says it is not well-defined..
The reason being that you could just do something like this:

ode:=diff(f(x),x) = -4*x/(x^2+1)^3 + 1/sqrt(f(x))-1-x^2

where the last 3 terms add up to zero for your given function.

If on the other hand the ode is required to be autonomous it seems that the question is more interesting.
If possible we want to find a function F of one variable so that

diff(y(x),x) = F(y(x))

has y(x)=f(x) as a solution at least on some interval I.

Take your example as an illustration.

restart;
f:=x->1/(1+x^2)^2;
solve(y=f(x),x) assuming y>0,y<=1;
finv:=y->sqrt(sqrt(y)*(1-sqrt(y)))/sqrt(y);
subs(x=finv(y),diff(f(x),x));
simplify(%) assuming y>0,y<=1;
F:=unapply(%,y);
F(f(x));
simplify(%) assuming x>=0;
ode:=diff(g(x),x)=F(g(x));
#dsolve({ode,g(0)=f(0)}); #Non-uniqueness at y=1=f(0), so choose another x0
dsolve({ode,g(1/2)=f(1/2)});
odetest(g(x)=f(x),ode) assuming x>=0;
#
## Just to illustrate the non-uniquess at y=1
DEtools[DEplot](ode,g(x),x=-2..5,[seq(g(1/2)=f(1/2)+.2*i,i=-3..2)],color=gray,linecolor=blue,thickness=1,stepsize=.05); p1:=%:
plot(f(x),x=-1..5,thickness=2); p2:=%:
plots:-display(p1,p2);

I changed your integration variable from x to t since it is confusing to have the integration variable and one of the limits bear the same name.

A:=-(5/8)*x*Int(f(t), t = 0 .. x)+(5/24)*Int(f(t), t = 1/3 .. x)-(1/4)*x*Int(t*f(t), t = 1/3 .. x);

combine(A);
evalindets(%,specfunc(anything,Int),IntegrationTools:-Split,1/3);
value(%); #To remove integrals that are clearly zero.
combine(%);
convert(%,list);
L:=map(`[]`@op,%);

g:=piecewise(t>=lhs(rhs(L[1,2])) and t<rhs(rhs(L[1,2])),L[1,1],t>=lhs(rhs(L[2,2])) and t<rhs(rhs(L[2,2])),L[2,1]);

I suppose that this could be automated, but I'm not going to do that.