@Rouben Rostamian  Sorry if i am being slow, but if the diagonalizable matrix we start with is real, then it can be expressed as PDP^{-1}.  In that case, the Ps are real as well, so the transpose of P is equivalent to the inverse.  A is real and symmetric so it can diagonalised with orthogonal matrices.

Or is is the case that P has to be unitary U in the decomposition for the case of the matrix exponential to work?  If that is the case I think my previous calculation is wrong, as I just showed that it was possible to decompose A into PDP^{-1}, where D is diagonal and P^{-1} is the transpose of P, and then used that to calculate the matrix exponential I needed, which looks similar to the one you have above.  Maybe this was just wishful thinking, as the transpose is obviously easier to calculate than the inverse.

@Rouben Rostamian  I'm not sure I follow, I have multiplied lambda with the transpose of the matrix U you have given, then multiplied by U, but I do not get A back.  Could you show explicitly how they multiply together to get A?

@vv I have the Fourier series, and I have another representation of cos (\alpha x).  This is the Fourier series for cos (\alpha x)/sin(\alpha \pi) and I used a different Fourier series to show that this gives another represention, which is superior to the approximation which the direct Fourier series gives.  However, I need to quantify the difference in error.  I hope that's clear.

@Rouben Rostamian   I am not sure if this is a full answer, since it is using the dsolve command and I'm not sure that the exercise wants you to do that.

To clarify, I proved for the Green's functon for the differential equation and specified conditions that G(x,u), diff(G(x),x) and diff(G(x),x$2) are continuous at x = u and also that in the limit as epsilon tends to 0, that diff(G(x),x$3) evaluated from x = u - epsilon to x = u + epsilon is equal to 1.

This was done using the delta function (since it's a fourth-order equation) and as you rightly say the related differential equation was to have the relevant Dirac delta function on the RHS.  As I say, I am pretty sure that using dsolve command and then simplifying is not what the exercise is asking for, as the exercises are fairly rigorous.

@Rouben Rostamian  Would you mind elaborating how this demonstrates the piecewise definiton for the Green function?  Sorry, I am a bit new to Green's functions, it seems like you are just assuming the definition which needs to be proved.

@Markiyan Hirnyk Could you elaborate on how this is showing the piecewise definition of G?

I was not familiar with Mr Baker's work (having not studied a lot of number theory) but I've looked it up and seems extremely interesting (and important), also that he did his degree where I am studying at the moment.  I hadn't actually realised that any Fields Medal-winning mathematicians had studied at my university, thanks for sharing the news of his passing.

@_Maxim_ Apologies, but could you explain the process a bit more explicitly as I still don't really understand what is going on and I am new to numerical analysis techniques, how is what you described helping me to calculate the coefficients c_k and P(n)?  Where does the b - h come from in the data points?

@_Maxim_ Thanks, I'll take a look again today.

@Torre Thanks a lot, I have done this now but seem to have got back a list which is much bigger than expected.  Some of these look to be duplicates, so there will be something for dz dz and then another value for dz dz.  Is everything I have obtained definitely a Christoffel symbol for the metric or could something else have gone wrong?  I have 30 values in total which just seems more than I was expecting, so wondering if it is normal to have some redundancy or repetition.

I have attached the updated code if you can maybe help me interpret what I am seeing:

FLRW_Metric.mw
 

Download FLRW_Metric.mw

@vv I have seen it stated in my Maple textbook that we can find an approximate value for S(n) in the limit by evaluating the second sum in your above equation in terms of the Riemann zeta function, is that what you are referring to with:

asympt(Sum(1/k^(11/10), k=1..n),n);

I have tried to input this sum into a procedure to observe the convergence but am getting error messages about the ; being unexpected.

S:= n -> evalf(Sum((1/k^0.1)*(sin(1/k)-1/k), k=1..n) + asympt(Sum(1/k^(11/10), k=1..n),n);));

Thanks for your reply, but I was explicitly looking for a procedure which would enable me to evaluate the sum, rather than just evaluating it by a direct command.

@Kitonum Can I also ask something else: one of my textbooks mentions this route as being part of the proof of the identity 

ln(sinh(x)^2 + cosh(x)^2 approx equal to (6*x^2(2*x^2 + 5))/(16*x^2 + 15)

Obviously once we have the values of the coefficients substituted in I can see how we can get to:

ln(sinh(x)^2 + cosh(x)^2 = (y(x)*(2*x^2 + 5))/(16*x^2 + 15)

but slightly stumped as to where the 6*x^2 comes from.

@acer I have managed to plot the graph I was asking about (although I did not really clarify what I had in mind) and obtained the perihelion values.  One last thing: I have looked up suitable parameter values and initial conditions appropriate for calculating the precession of Mercury.  How would I use these with your code to calculate the precession of the perihelion of Mercury?

@acer Thanks a lot, I'll take a look at that.  Will it be necessary to have the animations in order to determine precession of each revolution, or can I effectively have 'freeze-frames' after every 2*pi so that the perihelion can be determined each time?  Also given the solution r(phi) could you advise how I would use Maple to create a plot of the orbit in polar coordinates for r(phi) against r for a range of phi from 0 to 2*pi.