What is a "circular rod"?  I assume that you mean a rod with a circular cross-section.  That, however, is in conflict with your "1D" requirement.  If it is 1D, how does it matter that the cross section is circular, square, or whatever?

Anyway, if my interpretation of your question is correct, then you are looking at heat conduction in a line segment, let's say 0 < x < L.  The heat equation for the temperature u(x,t) is u_t = u_xx, where the x and t subscripts indicate derivatives with respect to the space x and the time t.  The initial condition is u(x,0)=u₀(x).  You need boundary conditions at x=0, and x=L.  Those depend on the details of what you want to solve.

Your statement "with a heat source at its base" is ambiguous.  Is the temperature of that heat source known?  If so, then the boundary condition at x=0 will be u(0,t)=A(t), where A(t) is the known temperature at x=0.  Alternatively, you may not know the temperature at x=0, but you may know the heat flux, that is, how much heat per unit time is injected into the rod. In that case the boundary condition would be u_x(0,t)=-B(t)., or in Maple notation, D[1](u)(0,t)=-B(t).  Similar considerations apply to the boundary at x=L, but you haven't said anything about that.

This may do what you want:

subs(tmp=1-theta, subs(theta=1-tmp, e1_1));

restart;
with(plottools): with(plots):
Obj := display(dodecahedron([0,0,0], 1));

# PQ is the axis of rotation
P := [1,-1,-1]: Q := [0,2,1/2]:
Axis := pointplot3d([P,Q], connect, color=red, thickness=5):

nframes := 30:
animate(rotate, [Obj, t, [P,Q]], t=0..2*Pi - 2*Pi/(nframes-1),
    frames=nframes, background=Axis,
    scaling=constrained, paraminfo=false);

restart;
de := diff(p(h),h)=A/(B+C*p(h));
ic := p(h0)=p1;
dsolve({de,ic}, p(h), implicit);

Thus, the problem has been reduced to solving a quadratic in the unknown p(h).  If you have some extra information about the problem's constants, then you may pick the desired solution from it.

Answer to question number 1: End each equation with a semicolon.

Answer to question number 2: Write the equations side-by-side, separate them with commas.

Do this:

sol := rsolve({f(n)=f(n-1)+10*f(n-2), f(1)=1, f(2)=2},f(k));
simplify([seq(sol, k=1..10)]);

The result is
[1, 2, 12, 32, 152, 472, 1992, 6712, 26632, 93752]

Solve the PDE numerically:
restart;
pde := diff(u(t,x),t)=diff(u(t,x),x,x);
ic := u(0,x)=sin(x);
bc := D[2](u)(t,-1)=0, D[2](u)(t,1)=0;
pdsol := pdsolve(pde, {ic, bc}, numeric);

Plot the solution at time t=0.2:
pdsol:-plot(t=0.2);

Animate the solution for a sequence of times:
plots:-animate(pdsol:-plot, [t=time], time=0..2,
   font=[Times,20], axesfont=[Default,12], axes=normal);

You should end every statement in Maple with a semicolon.  That includes compound statements with "if", as in:
if a > 1 then
  some statement;
  another statement;
end if;     # note the terminating semicolon of the "if" statement

In certain cases the semicolon may be optionally omitted but you don't have to.  Always terminate your statements with semicolons and you will be alright.

I can't tell what the inner workings of Maple may be on this, but there seems to be some good excuse to justify its behavior.

Maple has no trouble in solving that implicit solution for y(x):
solve(sol, y(x));

However, the structure of the solution is displayed much more clearly if we specify an initial condition:
sol := dsolve({ode, y(0)=1}, y(x));
plot(rhs~([sol]), x=0..4);

We see that the initial value problem is ill-posed—the initial condition does not determine the solution uniquely!  The y(0)=1 is not significant; any initial condition (other than y(0)=0) leads to non-unique solutions.  Thus, asking for a solution y(x) is not a meaningful request, and that's perhaps the reason why Maple stays on the cautious side.

This worksheet mw.mw shows how.  In it I have corrected a typo (a missing negative sign) in your calculations.

That said, it's not clear to me why you bother with Rayleigh-Ritz which was invented way back when, before the age of computers.  A modern numerical solver can obtain a more accurate solution, more easily.

I don't know what the azimuth and polar angles are.  In the attached worksheet I have used the longitude and colatitude angles which are normally used in spherical coordinates.  Perhaps they mean the same thing.

spherical-cap.mw

The problem with your plot command is that it does not refer to sol which contains the numerical solution of the differential equation. What you want is
   
or equivalently
   

You may add time range and color options, as needed.

One solves an equation.  One evaluates an expression.  I think you mean evaluate, not solve.

Anyway, here it is:

evalf(BesselI(0, 0.012*(-1)^(1/4)*sqrt(1000)));
        

I am guessing that you encoded ϕ_c as ϕ[c]. That's the source of your trouble. Instead, encode it as ϕ__c, that is, ϕ, followed by two underscores, followed by c.

Since you have not specified the parameter values, I will take your inequality to be:

and let q = 3.

restart;
eq := 1 + 2*q - 2*K*(1+q)/N;
q := 3;
tmp := plot3d(eq, K=0..1, N=0..1, view=0..8,
   axes=boxed, labels=[K, N, Gamma],
   filled=[style=surface, transparency=0.2]);
plots:-display(tmp, view=0.01..8);