I haven't worked on your problem simply because I don't like its input format. I know (from your emails) that you've read my extensive Post on boundary layer flow problems, "Numerically solving BVPs that have many parameters", and that the output (plots) that you want are very similar to my output in that Post. Please download the worksheet attached to that Post. Following my format, edit the beginning of that worksheet to use your ODEs, parameters, and boundary conditions. Do this using the Maple Input (1D input) like I have used. If you don't use that format, I won't read it.

Then enter into that worksheet the specifications for the output that you currently have in your attached .docx file. You can enter these into ordinary text fields. Then upload the edited worksheet to this thread.

If you need to discuss anything, please put it in this thread rather than sending me email.

@nm I wonder why your simplify doesn't cancel the 2 in the second term. Mine does.

@Carl Love My procedure above places the denominator of a numeric fraction coefficient in the numerator, which probably isn't wanted. Here's an updated procedure that handles them:

Apparent_Numer_Denom:= proc(e::`*`)
local F,R,N,D;
    (F,R):= selectremove(type, e, fraction); 
    (D,N):= selectremove(type, R, anything^negative);
    (numer(F)*N, denom(F)/D)
end proc:

Apparent_Numer_Denom(expr2); #Acer's expr2

@Anthrazit Applying evalb to the condition in an if statement is superfluous; it's done automatically.

@Test007 I didn't "come up with" any of the constants. The symbolic work is all done by dsolve, evalc, and limit.

The same process could be applied to the constant-coefficient ODE that you propose. Below, I treat the monic case, which covers all the cases with a second derivative:
 

Download ConstCoeffODE.mw

@Test007 I suppose that the feeling is that that situation can always be handled with eval. 

@tomleslie A procedure's behavior should not rely on the user having issued a certain with command at the top level, simply because the user may not have.

@nm Simpler:

add(ArrayTools:-Diagonal(A)^~2);

@Scot Gould As you proposed in an earlier version of the above Reply, the potential problem of using global t can be corrected with quotes:

drdt:= t1-> eval(diff(r('t'),'t'), 't'= t1):

Without the quotes, the following would give an error message, even if the assignment to t were totally unrelated to its use with r (such is always the risk with global variables):

t:= 3;
drdt(0);

@Scot Gould There are two problems with your derivative: 

• The symbolic derivative is recomputed for every invocation of drdt, including all the numeric points of the plot (if plotted in functional form).
• If global t is assigned a non-name value, the derivative won't work. 

@Carl Love Addendum: Using blocks= [u,v] should also work. It gives a different ranking of the derivatives. I think that the fact that it doesn't work should be considered a bug. Here is what the help page ?DifferentialRing says about the different rankings:

• Together, the lists blocks and derivations provide the ranking of R, that is, an ordering for the derivatives of the dependent variables with respect to the independent variables. In brief, the elements in blocks (dependent variables or sublists of them) are ranked between each other using an elimination ranking and those within a sublist are ranked used an orderly ranking.
  An elimination ranking is one where derivatives of leftmost dependent variables, or of all the dependent variables in a leftmost sublist within blocks, rank higher than derivatives of variables to their right. An orderly ranking is one where derivatives of higher differential order rank higher than derivatives of lower differential order, and if the order is the same, then derivatives of leftmost variables rank higher.

The precise bug that you encountered comes from a construction of the form a^b^c, which is always an error in Maple because it doesn't disambiguate between (a^b)^c and a^(b^c).

@janhardo By symmetry, the case A <> 0 and B=0 is equivalent to A=0 and B <> 0 with the roles of x and y switched. So, this case can be handled in your code by ans:= conics(0, A, D, C, E).

@nm Assumptions of equations specifying a value for a variable such as lambda=1 never work.

@emendes Just remove the one line that I specified. I don't think any other changes are needed.

Note that arg1 is not an argument to mainproc.

@janhardo The formula that you've used is sufficient for conics whose axes are parallel to the coordinate axes. For other conics, you need a term with both x and y. Conics include the degenerate cases of intersecting lines and a single point.