be able to solve it numerically, provided

1. you supply numeric values for the parameters q1, q2, q3, q4
2. you supply three initial/boundary conditions for the ODE

In the attached, I just use some random values for thes quantities

Download numODE.mw

You could start by typing

?How do I solve an ordinary differential equation

at the Maple prompt. If yoou need more info, check the Maple help for the command dsolve()

values for 'm' and 'l' (I just used 1 for both), then you could try DEplot3d() as shown in the attached. I just guessed a value for the range of the dependent variable p(t).

For some reason, this particular plot won't render on this site, by I can assure you that this code

VF.mw

  restart;
  with(DETools):
  sys := { diff(r(t),t)=p(t)/m,
           diff(p(t),t)=l^2/(m*r(t)^3)-n*k*r(t)^(n-1),
           diff(phi(t),t)=l/(m*r(t)^2)
         };
  sys1:=subs({n=1,k=1, m=1, l=1},sys);
  p1:=DEplot3d( sys1,
                [r(t), p(t), phi(t)],
                t=0..30,
                r=0..10,
                p=0..2*Pi,
                phi=0..Pi,
                arrows=cheap,
                dirfield=[5,5,5]
              );

will produce this plot

just your misunderstanding of how a plot command "works". When plotting some expression, the plot() command will evaluate the expression at ~200 points (by default), and then produce a "straight line" between adjacent points. The fact that you "see" this as a smooth "curve" is an optical illusion.

In your first example the points are only about 1.8 degrees apart (360/200) by default. In your second example, the points will be 3600000/200 degrees apart (by default) - so 18000degrees or about 50 complete revolutions. - so essentially you end up drawing a straight line from some point on the cicumference of a circle to some other (more or less random) point on the circumference. These lines will fill up the interior of the circle very quickly. You can get around this by increasing the option numpoints drastically and then wait a while for your plot to appear

@mmcdara 

the styles you want using conventional 'labels' and 'labelfont' options, see the attached

NB the double square brackets around the unit only appear on this site - they do not occur in the original worksheet

Download labs.mw

1. Two of your constraints are on the form s<p<q. Each of these has to be split into two simple inequalities, ie s<p, p<q. This fix is trivial
2. Both the expression you want to minimise and one of the constraints contain terms such as 91201.82720*(-0.05*i + 0.04)*sqrt(-i), - 912.0182718*i/sqrt(-i). Since 'i' is restricted to range 0..1, the sqrt(-i) factor in both of these terms is going to produce a complex value, which the optimizer cannot handle. ( All optimizers essentially work on the basis of determining whether an expression using the current variable values is greater than or less than the expression with the previous set of variables. However for complex values, the term "greater than or less than" has no meaning. So you have to ask, why do complex values occur?? Is this expected?
3. If the existence of complex values is unexpected then there must be an error in yur derivation somewhere. These complex valuesappear to arise from the evaluation of the quantities 'Q' and 'Qr'

I suggest you read the help page at ?pdsolve/numeric, very carefully

The attached at least executes, although I can't say I'm wildly convinced by the answer which Maple provides

Inline displlay of worksheets on this site appears to be broken - again!

Download pde.mw.

the command line (which *seems* to be true), then you  can add

 lprint(diff(ln(x),x));

to the end of your probem.mpl file. This will produce 1/x

If you don't need/want to show the steps, then you can just change the problem.mpl file to

  diff(ln(x),x);

which will also produce 1/x;

Your integral can be solved without resdtriction(s)- see the attached

Download simpleInt.mw

as in the attached ( a "toy" example)

Download nyplt.mw

to set the x and y-ranges appropriately - yours are much too small.

See the attached

Download ellip.mw

hack something in the Finance() package to produce what you want, but since there is a simple closed form expressionj for a standard mortgage (see about halfway down on  https://en.wikipedia.org/wiki/Mortgage_loan), it is pretty trivial to do in an explicit manner.

In the attached, there are two methods of doing this calculation, The first uses a recursive function call, basically keep adding interest and subtracting payments, until the duration is up, then solve for the payment which makes the outstanding loan zero at the end of the duration. The second just uses the formula from the above Wikipedia page - you can google for the derivation

You can give students extra points if they manage to get the output to appear (as in the attached), by using the  settings under the menu Format->NumericFormatting

Download mort.mw

is to use isolate() as in the attached

Download isol.mw

the ratio of circumference to diameter is Pi not pi - so the attached will work

Download plotA.mw

completely clear to me what you want since you state that

i need the solution  for Y=0 and Y=1

but you also seem to require the integral expression

(int(U(Y)*Theta(Y), Y = 0 .. 1))/(int(U(Y), Y = 0 .. 1));

The attached covers both of these requirements.

Download anODE.mw