Of course the more specific you are about the function h(x) and the initial/boundary condition, the more detailed the solution will be. The help page is found online here.

Download Bernoulli.mw

Here's my take on the derivatives, since you explicitly asked about D

Download Derivatives.mw

For approximations, you always trade off compactness for accuracy, so I don't see why you don't just plot derivatives. To solve any confusion about which Root is used, you can check with fdiff, which does it by numerical function evaluations of the RootOf. In all the cases here, it is the same as a derivative using D

MaPal93approx.mw

For your directed acyclic graph, the algorithm is relatively simple and can be be done with a compiled Maple routine that runs in a reasonable time (3 min on my old laptop).

Download BacktrackDAGPaths.mw

DAGPaths:=proc(G::Graph,	# directed acyclic graph - self loops and duplicate edges (if multigraph) are ignored.
			V1, V2,	# start and end vertices to find paths between.
			{minlength::nonnegint := 0, 	# (optional) minimum length of path (measured in number of edges).
			maxlength::nonnegint := -1})	# (optional) maximum length of path - default is length of longestpath.
			# outputs an npaths x (maxlength+1) Matrix (datatype = integer[4]) that contains the paths.
			# entries are positive integers corresponding to the vertices.
			# use infolevel[DAGPaths] := 2 to output additional information. 
	description "finds all paths between two vertices in a directed acyclic graph";
	uses GraphTheory;
	local P, S, A, Adj, n, i, v1, v2, GG, G2, VerticesG, maxpathlength, minlengthi::integer[4], maxlengthi::integer[4],
		maxpaths, maxdeg, pathptr, path, paths, npaths, N, j, src::integer[4], dest::integer[4], newvertices;	
	# check arguments and set up 
	if not IsDirected(G) or not IsAcyclic(G) then error "graph must be directed and acyclic" end if;
	VerticesG := Vertices(G);
	if not member(V1,VerticesG,'v1') then error "Vertex %1 not in graph", V1 end if;
	if V1 = V2 and minlength = 0 then
		return Matrix(1, 1, [1], 'datatype' = integer[4])
	elif V1 = V2 and minlength > 0 then
		userinfo(2, procname, "no paths between %1 and %2 with length greater than %3", V1, V2, minlength); 
		return Matrix(0, 1, 'datatype' = integer[4])
	end if;
	if not member(V2, VerticesG, 'v2') then error "Vertex %1 not in graph", V2 end if;
	if OutDegree(G, V1) = 0 then error "start vertex has no departures" end if;
	if InDegree(G, V2) = 0 then error "end vertex has no arrivals" end if;
	if HasSelfLoop(G) or IsMultigraph(G) then
		userinfo(2, procname, "ignoring duplicate edges and self-loops")
	end if;	
	n := nops(VerticesG);
	GG := UnderlyingGraph(G, 'directed' = true, 'selfloops' = false, 'multigraph' = false);
  	# Delete departures from last vertex
  	if OutDegree(GG, V2) = 1 then
  		DeleteArc(GG, [V2, Departures(G, V2)[]])
  	elif OutDegree(GG, V2) > 1 then
  		DeleteArc(GG, [seq([V2, i], i in Departures(G, V2))])
  	end if;
  	# and delete arrivals to first vertex
  	if InDegree(GG, V1) = 1 then
  		DeleteArc(GG, [V1, Arrivals(G ,V1)[]])
  	elif InDegree(GG, V1) > 1 then
  		DeleteArc(GG, [seq([V1, i], i in Arrivals(G, V1))])
  	end if;
  	# work only on component that V1 is in for path length calc.
  	# (could be ommitted but calc is faster with smaller adjacency matrix)
  	newvertices := select(has, ConnectedComponents(GG), V1)[];
  	if not member(V2, newvertices) then
  		userinfo(2, procname, "no paths between %1 and %2", V1, V2); 
		return Matrix(0, 1, 'datatype' = integer[4])
	end if; 
	# Find length of longest path
	G2 := InducedSubgraph(GG, newvertices);
	Adj := AdjacencyMatrix(G2);
	maxpathlength := nops(BellmanFordAlgorithm(Graph(-Adj), v1, v2)[1]) - 1;
	userinfo(2, procname, "longest path length is %1", maxpathlength);
	if maxlength = -1 then
		maxlengthi := maxpathlength
	else
		maxlengthi := min(maxlength,maxpathlength);
	end if;
	if minlength > maxlengthi then
		userinfo(2, procname, "no paths with minlength %1 greater than maxlength %2", minlength, maxlengthi);
		return Matrix(0, 1, 'datatype' = integer[4])
	else
		minlengthi := minlength
	end if;
	userinfo(2, procname, "calculating paths with length between %1 and %2", minlengthi, maxlengthi);	
	# calculate number of paths from adjacency matrix since # walks = # paths
	P := Adj[v1];
	for i from 2 to minlength do
	P := P.Adj;
	end do:
	S := P:
	for i from minlength + 1 to maxlengthi do
  		P := P.Adj;
  		S := S + P;
	end do:
	maxpaths := S[v2];
	userinfo(2, procname, "calculating %1 paths... ", maxpaths);
	A := op(4,GG):
	# Construct an n x maxdeg Matrix equivalent to op(4,GG) 
	maxdeg := max(map(nops, A));
	N := Matrix(n, maxdeg, 'datatype' = integer[4]):
	for i to n do
		for j to nops(A[i]) do
			N[i,j]:=op(j,A[i])
		end do;
	end do;
	# allocate space and initialize
	src := v1;
	dest := v2;
	path := Vector(maxpathlength + 1, {1 = src}, 'datatype' = integer[4]); # +1 from edges to verts
	paths := Matrix(maxpaths, maxlengthi + 1, 'datatype' = integer[4]);
	npaths := Vector([0], 'datatype' = integer[4]);
	pathptr := 1;
	# call DFS recursively.
	DFS(src, dest, npaths, pathptr, path, paths, minlengthi, maxlengthi, N, maxdeg);
	if npaths[1] = 0 then
		userinfo(2, procname, "no paths between %1 and %2 with lengths between %3 and %4", V1, V2, minlengthi, maxlengthi)
	elif npaths[1] = maxpaths then
		userinfo(2, procname, "found expected number of paths")
	else
		userinfo(2, procname, "found %1 paths - unexpected result", npaths[1])
	end if; 
	# output results
	paths;
end proc:
#
# Compilable recursive Depth-First Search routine - compiled below
#
	DFSM:=proc(src::integer[4], dest::integer[4],
			npaths::Vector(datatype = integer[4]), pathptr::integer[4], path::Vector(datatype = integer[4]),
			paths::Matrix(datatype = integer[4]), minlength::integer[4], maxlength::integer[4],
			N::Matrix(datatype = integer[4]), maxdeg::integer[4])
		local v::integer[4], j::integer[4];
		if src = dest then
  			if pathptr >= minlength + 1 and pathptr <= maxlength + 1 then
    				npaths[1] := npaths[1] + 1;
	    			for j to pathptr do	
					paths[npaths[1], j] := path[j];
				end do;
			end if
		else
			for j to maxdeg do
				v:= N[src, j];
				if v = 0 then break end if;     
				path[pathptr + 1] := v;
				procname(v, dest, npaths, pathptr + 1, path, paths, minlength, maxlength, N, maxdeg);
			end do;
		end if;
		NULL;
	end proc;
#
DFS := Compiler:-Compile(DFSM);

Download infix.mw

(I didn't understand what you meant by simplification if not the above.)

This is fine for simple units but is fairly limited for compound units which are just formatted as Maple expressions with the "*" etc. And for printing formatted strings, symbols such as Angstroms are expressed as their html code, so it would take a lot of work to format unit strings for the general case.

There are lines in the worksheet that are not displayed in Mapleprmes, e.g.,

printf("The voltage is %a %a", NumberUnit(x))

Download PrintUnits.mw

Here's a second version you may prefer

PrintUnits2.mw

printf("The electric field in base units is %5.2f %s", NumberUnit(x))

The electric field in base units is  5.00 m kg / A s^3

The plot option useunits means the units are interpreted correctly. (And your f2 in 1D was missing a multiplication.)

FunctionUnits.mw

It could probably be made more efficient but I wasn't sure I understood the descriptions and all the assumptions.

Download sort2.mw

seq(NumberTheory:-CalkinWilfSequence(i),i=1..10)

1, 1/2, 2, 1/3, 3/2, 2/3, 3, 1/4, 4/3, 3/5

Here's one way.

Download TabPropQ.mw

(You may want to exclude 1 from k.)

This may be a helpful, though oversimplified example. As for visualization, once you have a solution, each lambda can be plotted in a 3-D plot against two of the parameters for fixed values of the other parameters. It may be possible by some sort of "dimensional" analysis to find combinations of the variables that reduce the number of variables/parameters.

Download symmetry.mw

Heres another approach for more recent Maple versions in 1-D only. It could be worked up into a procedure as @mmcdara did, but I am short of time now. I think it should be reasonably efficient but I didn't do any tests.

Download search.mw

As @Carl Love says, it is common to get these small imaginary parts. They can be removed with fnormal and simplify(...,zero). Your case is unusual in that the digits for fnormal needs to be significantly different from the digits for the calculation - I don't recall ever having to use the digits option to fnormal before. Presumably the inaccuracies in the dilog, etc functions accumulate errors.

[Edit: the imaginary part looks suspiciously like Pi times a power of 10, so perhaps there is more going on here.]

Download fnormal.mw

for m from 0 to 5 do
  printf("% 5d\n",<seq(binomial(m,i),i=0..m)>);
end do;

As in @nm's answer, I'm not quite there yet, but have resolved the issue of which is the correct root

Download RootOf.mw

You can work around this by making a new context for the are unit in which any SI prefix is allowed.

Download Areas.mw