This website still fails to display uploaded worksheets.  Is anyone looking after this?

I will write L for the line that you have called the axis of the bundle.

Let a be a (nonzero) vector along L, and let p the position vector of a given point on L.

Pick an arbitrary nonzero vector s which iis not parallel to a.  Any such vector will do.  Then the vector b = a x s would be orthogonal to a, and the vector c = a x b would be orthogonal to both a and b.

Normalize the vectors a, b, and c, that is, divide them by their own lengths.  Then for any angle θ, the unit vector d = b cos θ + c sin θ is in the plane of b and c, and therefore perpendicular to a. It follows that a parametric equation for the planes spanned by a and d and going through the point p is
r = p + u a + v d.
where u and v are the parameters of the plane.  As θ varies, you obtain the bundle of planes that you are looking for.

See the attached worksheet for a worked-out example.  The following illustration is extracted from that worksheet.  The point p and the vector a are shown in black.

Download mw2.mw

Here is a solution.  I have taken advantage of Maple's built-in Euler-Lagrange equations to derive the geodesic differential equation from scratch.  Alternatively, I could have just typed in the differential equation from a textbook.

Unfortunately this site is broken (again!) and won't allow displaying Maple worksheets.  Click on the link to download it. This figure is extracted from the worksheet.

Download geodesic-on-torus.mw

@MAXR, You are using Maple's default "Document Mode" for input which is not easy to work with.  Because of that, your document is full of syntax errors.  As sand15 suggested, you will have a better chance of success if you switch from "Document Mode" to "Worksheet Mode".  Follow the instructions here to swith your Maple 17 to "Worksheet Mode".  Myself never use the "Document Mode".

Your document is so full of syntax errors that it's not worth the effort to rescue it. You will be better off by starting with a  blank worksheet and re-entering your code.

I have started it for you in the attached worksheet.  Start from there and enter the rest.    Observe that the "Worksheet Mode" uses what is called a "1 D" method of entering math.  You should be able to figure out how by looking at that sample.

Also note that as I had pointed out in a previous reply that you cannot use both lambda and lambda[1] in your code because they conflict.  Change lambda to lambda[0] everywhere in your code as you type it.

And an advice:  Rather than entering your equations all at once, enter them one at a time.  Again, follow the example in the worksheet.

In the boundary condition, I have changed your diff(theta(0), eta) to D(theta)(0).  In your version it evaluates to zero because theta(0) is just a constant, so its derivative is zero, which is not what you want.

Note that I have broken long equations lines into multiple lines for the sake of readability.  That's not essential but makes it easier to spot errors. To break a line at the cursor's position, press "Shift Enter".

I have followed your method of entering the third derivative as diff(f(eta), eta, eta, eta).  I would enter that as diff(f(eta), eta$3).  Use whichever form you prefer.

Even after fixing all the syntax errors, there is no guarantee that Maple will find a solution.  A boundary value problem may have no solutions at all unless the data is specified just right.

Note: Be sure to verify that what I have entered in this worksheet is correct. 

Download maple-worksheet.mw

The answers presented thus far do correctly show that for any n within a finite range, such as n=0..25, there is an integer m that satisfies the requirements.  But that doesn't prove that the conclusion holds for arbitrary n. The following worksheet presents a proof.

Download mw.mw

It works for me.  According the formula that you have provided, D(ln@cosh) should evaluate to

D(ln) o cosh 
times
D(cosh)

But D(ln) evaluates to z->1/z  and D(cosh) evaluates to sinh.  Therefore Maple duly returns

((z -> 1/z)@cosh) sinh

which is the correct answer.  

Earl, here is how I would solve the problem.

We have a curve (a track) in the vertical plane, and a bead that slides along the track.  The coefficient of friction is μ.  We wish to describe the bead's motion.

It's possible to derive the equation of motion entirely in Maple as you are attempting to do, however the reasoning is more transparent when done by the proverbial "paper and pencil" as in the attached sliding-with-friction.pdf.  The calculation is generic, in the sense that it applies to any track shape and any choice of parameters. 

In the attached Maple worksheet I have implemented the resulting equation, and then illustrated by applying it to the specific case of a circular track.  The bead slides down the circle, goes past the bottom, and oscillates a few times before coming to a standstill abruptly in finite time due to friction.  You may change the shape of the track and other parameters as you wish.  They are all defined in the statement params := ... in the worksheet.

>  restart; >  with(plots): >  Typesetting:-Settings(typesetdot=true): >  r := s -> < x(s), y(s) >; (1) The tangent vector >  T := diff(r(s),s); (2) The length of the tangent vector >  T_len := sqrt(T^+ . T); (3) The unit tangent >  T / T_len: convert(%, D); tau := unapply(%, s): (4) The unit normal >  <-tau(s)[2], tau(s)[1]>: convert(%, D); nu:= unapply(%, s): (5) The expression : >  diff(r(s),s): convert(%, D): rp := subs(s=s(t  ), %); (6) The expression : >  diff(r(s),s,s): convert(%, D): rpp := subs(s=s(t), %); (7) The length of rp >  rp_len := sqrt(rp^+ . rp); (8) The velocity >  v := diff(r(s(t)),t); (9) The gravitational pull: >  gvec := <0, -g>; (10) Here is Q (see documentation) >  simplify(tau(s(t))^+ . v); Q := mu*signum(%); (11) Newton's Law of Motion:   Mass times acceleration = applied forces See the documentation for the derivation. >  rp_len*diff(s(t),t,t)         + rpp^+ . ( tau(s(t)) + Q*nu(s(t)) ) * diff(s(t),t)^2         - gvec^+ . ( tau(s(t)) + Q*nu(s(t)) ): de := simplify(%); (12) >    Example   Everything up to this point is generic, as we have not made any assumptions on the shape of the curve or other parameters.  Here we work out a specific example. >  params := g=1, mu=1/10,         x = (s->  -cos(s) + 1),         y = (s->  -sin(s) + 1) ; (1.1) Here is what the track and the starting point look like: >  subs(params, [x(s), y(s), s=0..2*Pi]): track := display(         plot(%, color="Blue", thickness=3),         pointplot(eval([x(0),y(0)], {params}),                 symbol=solidbox, symbolsize=30,                 color="Green"), scaling=constrained); and here is the differential equation for this particular instance: >  subs(params, de): DE := simplify(%); (1.2) Solve the DE: >  dsol := dsolve({DE, s(0)=0, D(s)(0)=0},         numeric, output=operator); (1.3) Plot s(t): >  odeplot(dsol, [t,s(t)], t=0..20); Warning, cannot evaluate the solution further right of 16.400954, maxfun limit exceeded (see ?dsolve,maxfun for details)   The solution extends up to the time where the bob stops moving. We set tmax according to Maple's warning message above; >  tmax := 16.400954; (1.4) This proc displays the bob at a time t: >  frame := proc(t)         local S;         S := subs(dsol, s);         subs(params, [x(S(t)), y(S(t))]);         pointplot(%,                 symbol=solidcircle, symbolsize=30, color="Red"); end proc: Animate the bob's motion.  Friction stops the bob after a few oscillations. >  animate(frame, [t], t=0..tmax, background=track, scaling=constrained, frames=60); >

Download sliding-with-friction.mw

Maple's export to PDF has been broken for a long time, and despite my occasional prodding, it remains useless.

Try, for instance:

restart;
with(plottools):
plots:-display(disk([0, 0], 1), disk([1, 0], 1), color = [red, green], transparency = 0.8, scaling = constrained);

Right-lick on the image and select Export -> PDF, and save.

Viewing the saved result we see two glaring shortcomings:

1. The transparency is lost.  That's despite the fact that the image is encoded as a PDF 1.4 which has native support for transparency.
2. The graphics is embedded in a 612 x 792 pts (US letter size) "paper".  For the life of me, I have no idea how graphics embedded in a letter-size page can be of use to anyone.  Shouldn't the graphics be saved in dimensions compatible with its bounding box?  I seem to recall that Maple 11 used to do that properly.

As an alternative to PDF, you may want to try saving your graphics in the Encapsulated PostScript format.  That provides faithful vector graphics, and respects the image's bounding box.  It does not handle transparency because PostScript doesn't.

By the way, the export menu spells "Postscript" with lowercase "S".  The proper spelling should be with an uppercase S as that is a trademark.  That should be corrected.

For whatever it's worth, ChatCPT derives a very simple answer through elementary trigonometry.  I don't know how to show ChatCPT's derivation here, so I will just show the final result:

f := (n,x) -> n*x + 2*sum((n-k)/k*sin(k*x), k=1..n-1);

In case you are interested, I obtained that by asking ChatCPT:

integrate 1-cos(n*x))/(1-cos(x)) where n is an integer

@C_R The following worksheet outlines the technique for solving the initial boundary value problem corresponding to the boundary condition BC3 = u(a,z,t)=0.  The insulated boundary condition may be treated similarly.

The implementation of the technique is tedious but straightforward if you are familiar with separation of variables and eigenfunction expansion of solutions of PDEs.  janhardo has provided AI-generated details of some concrete cases.

Download heat-pde-outline.mw   <-- contains typos,  don't use!

Edit: 

I spotted a view careless typos in the originally posted worksheet.  These have been corrected in this version 2.

Download heat-pde-outline-ver2.mw

A thorough answer will depend on what you intend to do with the result.  Based on what you have shown so far, the following should serve your needs.

Download mw.mw

Edit:  Alternatively, you may simply do

restart;
e := Matrix(3,3, (i,j) -> 1/2*(diff(u[i](x[1],x[2],x[3]),x[j])
                               + diff(u[j](x[1],x[2],x[3]),x[i])));

To do what you want, specify the file name as ./data.xlsx, as in

restart;
with(ExcelTools):
Import("./data.xlsx");

This was tested and works as expected in Linux.  It should work in other operating systems as well.

You may solve a set of ODEs interfacing over any number of adjacent intervals through the method of shooting.  The idea is that you have the value of the solution at the left end, you leave the value of the slope there as an arbitrary parameter, and then apply fsolve() to determine the value of the parameter so that the solution hits a desired value at the right end.

The following worksheet illustrates how to solve such a problem over two adjacent intervals.  It may be extended in the obvious way to solve the interface problem over any number of adjacent intervals.

See if you can adapt this to solve your problem.

Download shooting.mw

ithprime(b) $b=1..19;

After a bit of algebra it is possible to reduce the original integral to the following but I don't know how to evaluate the result symbolically.  If anyone can figure out how, then I will post my calculations that lead to that result, otherwise it's not worth the trouble.

Download mw.mw

Here an implementation of the equations of motion derived in
https://www.youtube.com/watch?v=1h4zlDRjV0c
but you may want to have a look at the older Maple demos that the other posters have noted.  They are more illuminating.