@Kitonum Thanks for the answer.  I somehow did not think of applying expand().

This wroksheet demonstrates how to compute the motion of a ball that bounces on a hilly terrain:

The worksheet: bouncing-ball-on-hilly-terrain.mw

This generalizes my calculation from an earlier post where the ball was bouncing vertically on a flat ground.  In that calculation the gravity below the ground was set to zero in order to prevent the ball penetrating the ground.

Unfortunately that trick does not work in the non-flat case.  Actually, I don't quite understand why it worked at all in the flat case. I must have been wearing my Aristotelian Physics hat in believing that if there is no force, then the motion would cease.  But Newton has taught us that an object in motion remains in motion even when there is no force acting upon it.  This is amply demonstrated by setting the gravity to zero under the hilly ground and observing that it does not stop the ball from penetrating.

To get around that, I modified the model by letting the gravity pull down as it normally does above the ground, but to push up much more strongly from below the ground. That prevents the penetration.  The animation above show a sample result.  You may play around by changing the parameters in the worksheet.

@vv That's nice.  I knew how to calculate the individual arcs, but it had not occurred to me that it would be possible to put all the pieces together in a single formula.

@Preben Alsholm I calculated the stopping time and I obtained 3*sqrt(2), agreeing with yours.  In general, the stopping time corresponding to the initial conditions y(0)=1, D(y)(0)=0 is

where g is the acceleration due to gravity, and r is the speed reduction factor at an impact. As you have noted, there are infinitely many bounces during this finite time interval.  Some people call that the Zeno phenomenon.

OK, guys, inspired by your answers, I made yet another attempt for a solution, and this one works very nicely.

What we do is to let the acceleration of gravity, g, be a nonzero constant above the ground, and be zero under the ground.  Then everything works perfectly and the ball comes to a standstill on its own after about t=4 in agreement with Preben's observation.  This approach does not require the introduction of artificial "small threshold" numbers.

restart;
g := y -> piecewise(y>=0, 1, 0);
de := diff(y(t),t,t) = -g(y(t));
ic := y(0) = 1, (D(y))(0) = 0;
Events := [y(t), diff(y(t), t) = -.5*(diff(y(t), t))];
dsol := dsolve({de, ic}, numeric,
            events = [Events], range = 0 .. 6);
plots[odeplot](dsol, thickness = 3, color = red);
plots[odeplot](dsol, [t,diff(y(t),t)], thickness = 3, color = red);

@Preben Alsholm I agree that the help page is not very helpful.  Not only that, it is essentially impossible to reach unless you know exactly what you are looking for.  I have made a little note to myself about how to reach help regarding events, and I refer to it when I need help.  For everyone's benefit, here it is:

In Maple's search box enter
    Events for dsolve[numeric]    (or equivalently, at Maple's prompt enter ?dsolve,events)

For more examples, in Maple's search box enter
    How Do I Solve an Ordinary Differential Equation?

@Mariusz Iwaniuk Thanks.  This is a good idea. In effect, it says when the ball digs itself underground, then just keep it there.  Thus,  beyond the halt event, y(t) takes on a constant small negative value.

I had thought of this solution but did not try it because I was under the impression that once a halt event occurs, the computation stops and the solution is undefined beyond that.  Looking at your solution however, it appears that that's not the case.  Plotting the solution shows that the solution is well-defined beyond the halt point.  This may be in the documentation, by as Preben has noted, the documentation is not the easiest thing to read.

@acer Thanks.  This is a reasonably good solution.  In effect, it turns off the amplitude reduction mechanism once the amplitude falls below a certain threshold.  Beyond that the ball continues bouncing up and down forever with that small constant amplitude.

What you are saying here is either wrong, or the information provided is incomplete.

Here is my understanding of your question.

You have the equation

    a(x,y,t,u(x,y,t),...) = the_right_hand_side;

You expect the derivative of that equation with respect to x and to be:

    a[1](x,y,t,u(x,y,t),...) = the x derivative of the_right_hand_side.

That is certainly not true, unless you have some hidden assumptions which you have not revealed.

The derivative of the equation with respect to x will involve not only a[1], by also a[4], etc.

@weidade37211 It is difficult to tell what is wrong with your code without seeing the code.

@weidade37211 There is no requirement of a "closed form" f(x,y).  All we need is something that receives x and y and produces a value.  Its internals may be as complicated as you wish, and may contain any number of calls to fsolve().

@Jjjones98 I cannot tell what you are calculating here.  There are several functions involved in this discussion -- z(x), f(x), y(x), G(x,u).  Which of these are you calculating?  Or perhaps you are calculating something entirely different?

@Jjjones98 In my calculations the function f(x) was selected to satisfy the inhomogeneous boundary conditions only.  There was no requirement at all that f(x) were to satisfy the differential equation, although it turned out that it did, merely by accident, not by design.  In fact, in a more complicated case, for instance when the differential equation has non-constant coefficients, finding an f(x) that satisfies both the boundary conditions and the differential equation can be quite impractical.

Your formulation where you require f(x) to be "the solution of the homogeneous equation with inhomogeneous BCs" is unnecessarily restrictive although in the case of this essentially trivial differential equation does not pose a problem.

This question was asked and answered here:

https://www.mapleprimes.com/questions/224240-Perturbation-Theorie#comment247231

@Jjjones98 I will give you a sketch of the approach.  You should be able to fill in the details.  Ask again if you run into problems.

Your Green's function matches the homogeneous boundary conditions
y(0)=0, y'(0)=0, y(1)=0, y'(1)=0,
while the problem you wish to solve calls for the inhomogeneous boundary conditions
y(0)=1, y'(0)=1, y(1)=2, y'(1)=3.

Your first step is to reconcile these.

1. First, find a function, any function, that satisfies the inhomogeneous boundary conditions.  We have four conditions to satisfy, therefore the cubic function
   f := x -> a*x^3 + b*x^2 + c*x + d
   which has four undetermined coefficients is a good candidate.  A straightforward calculation shows that
   f := x -> 2*x^3 - 2*x^2 + x + 1.

2. Let's look at the boundary value problem that you wish to solve:
   y''''(x) = sqrt(x),   y(0)=1, y'(0)=1, y(1)=2, y'(1)=3.
   Define z(x) = y(x) - f(x).  Since y(x) and f(x) both satisfy the inhomogeneous boundary conditions, their difference z(x) satisfies the homogeneous boundary conditions.  That's good news.

3. From z(x) = y(x) - f(x) we get y(x) = z(x) + f(x).  Plugging this into the differential equation we obtain z''''(x) + f''''(x) = sqrt(x).  But f(x) is a cubic, therefore its fourth derivative is zero.  We conclude that:
   z''''(x) = sqrt(x),    z(0)=0, z'(0)=0, z(1)=0, z'(1)=0.

4. We solve for z(x) by applying your Green's function.  You may do that in Maple, as in
   z(x) = int(G(x,u)*sqrt(u), u=0..1) assuming x > 0, x < 1;
   or by hand, if you are patient enough.  In either case you will obtain:
   z(x) = 16*x^(9/2)*(1/945)-8*x^3*(1/189)+8*x^2*(1/315).

5. Now we put z(x) and f(x) together to arrive at
   y(x) = z(x) + f(x) = -622*x^2*(1/315)+16*x^(9/2)*(1/945)+370*x^3*(1/189)+x+1
   which agrees with what you have posted.