I am trying to use Maple to compute the numerical solution to the Newton/Einstein orbital equations and then to plot r(phi) against phi in polar coordinates (obviously without the relativistic correction this should be a flawless ellipse).

I have seen that problems trying to solve ODEs/PDEs are fairly common with people who are new to Maple.  The problem I have is that when I define the ODE with the second derivative, it evaluates the derivative so when I try to use dsolve I get the error message that it was 'expecting' ode to be a differential equation.

ode := diff(1/r(phi), [phi$2]) + 1/r(phi)=GM/h^2;     [(where constants are defined such that GM/h^2 = 1)]
dr :=diff(r(phi),phi);

ics := r(0)=2/3, dr(0)=0;

dsolve([sys,ics]);

Also, I am not sure how to define one of the conditions when it is a derivative such that it can be included with sys, although I hade made an attempt.  The other examples I have seen just have functions defined at a particular value, rather than a derivative of the function.  When I try to use dsolve (regardless of the other problem), I get the message 'Error, (in dsolve) not a system with respect to the unknowns {r(0), (diff(r(phi), phi))(0)}'.  Once I have the solution I would like to plot r(phi) against phi for 0 < phi < 2pi to observe the full orbit and then do the same for 0 < phi < 10pi once I have repeated the numerical solution for the orbital equation with relativistic correction so that I can observe that the precession of the perihelion between each revolution is approximately constant. 

I am also not sure how to find the angular positions of the four perihelia after the first revolution, is there a command that will enable me to determine the distance from the origin of the orbit to the closest point on the orbital path for all of the revolutions over the range 0 < phi < 10pi, as I will then be able to make a comparison with each one to check that the precession is approximately constant.

I have a function of the form:

y:= x-> ((1 + ax + bx^2)/(1 + cx + dx^2))*(ln(sinh(x)^2 + cosh(x)^2)

I would like to know how I could use Maple to calculate the values that the constants a,b,c and d should take such that the expansion of the above function does not include powers of the order x^3, x^4, x^5 or x^6 ie. such that the powers are quadratic at most.

The trigonometric terms just to clarify are the square terms ie. sinh(x) * sinh(x), but that is how I have written them before with Maple.  Not sure if I have written it out correctly, but it is a fraction with the constants multiplying the natural log function whose argument is the sum of the squares of the hyperbolic trigonometric functions.

I am working on a simple PDE and am trying to find the values of the constant coefficients in the PDE such that a function of a particular form satisfies the equation.  The function is a sum of general functions of x, y and z in the argument: for example, one of them is $f(x-2y+z)$.  I then have a PDE with constant coefficients which contains a sum of some mixed and some single second partial derivatives of that function with constant coefficients in front of some of the derivatives.  Is there any relatively easy package with Maple which would enable me to find the values of those coefficients such that the function of the form I have will satisfy the PDE?  I can provide more exact details of the PDE if necessary but didn't want to get bogged down in details as interested in the principle.

I am attempting to use Maple to solve a matrix equation of the form aB^2+bB+cI=B^(−2), where B is a 3×3 matrix whose entries are constants which have been specified, and I is the 3×3 identity matrix (i.e. to find the values a, b and c which satisfy the equation).

My idea was to use a piece of code like

A:=Matrix([[7,4,-2],[4,7,5],[2,-3,8]]);
d:=Vector([[8],[5],[2]]); 

for matrix inversion to solve a system of three linear equations but that would assume that the right-hand side of the equation is a column vector when it is actually another 3×33×3 matrix. One other approach I took was to equate the left and right hand sides so that you add the three matrices on the LHS and then equate the resulting each entry with the corresponding entry of the matrix on the right hand side, however this gives me a set of 9 linear equations for 3 unknowns.  When I try to solve these with fsolve it says that there are more equations than there are unknowns.