Consider the difference between

simplify(subs(z = -z, erf(z)));
simplify(subs[eval](z = -z, erf(z)));

Help page for the 'subs' command does state

The action of substitution is not followed by evaluation by default

Step-by-step solution information is only available for limited number of commands within some 'Student' sub-packages.

None of these cover the solution of differential equations. Hence no step-by-step illustration is available for your requirement

dsolve() is used to solve differential equations. You appear to be producing a pair of (non-differential) equations involving C1 and Nu, which can be solved for all specified values of Br, by using fsolve(), as in the attached

doSols.mw

It takes a bit of juggling of ranges for implicitplot3d() calls, and 'view' regions for the final plot to get everything in one plot in a more or less visible manner. See the attached

Download plotprob.mw

Generally agree with Kitonum that re-executing OP's worksheet does not produce plots with a logarithmic x-axis. I asssume that OP has been using all sorts of "interactive" plotting commands from context menus to actually produce the displayed figures - none of whihc are visible, making dignosis a little more difficult. However OP's problem can be illustrated programmatically by using using the 'mode=log' option for 'axis[1]', in a simple odeplot() command. One can appreciate the same effect on a linear x-axis by setting the plotting range to 1e-4..10and the 'view' option to (say) [1e-04..2e-01, 0..1]

The problem arises because the adaptive plotting routines are selecting the points-to-plot 'linearly', but the mode=log option then forces the x-axis to display them 'logarithmically'. (I'm using the terms 'linearly and 'logarthmically' here rather inaccurately, just to explain the problem). The net effect of this process is that far too few points are generated with very small values of the independent variable 't', when the t-axis is subsequently set to display logarithmically. Hence one gets a sort of 'piecewise plot' effect for small values of the independent variable.

As the attached worksheet shows, there are (at least) two ways to fix this.

1. The first one is just to increase the value for ' numpoints' option in the odeplot() command - ie generate a lot more points. This means that the spacing between plotted points  for small values of the independent variable will be greatly reduced, and the 'piecewise' effect is correspondingly diminished. The drawback is that one one is now generating "too many" points for larger values of the independent variable, slowing down the plotting process. (NB this takes about 10secs on my machine, so just wait).
2. One *might* think that a *similar* effect could be obtained by setting the 'maxstep' option in the dsolve() command to something very small. However it *seems* that the plotting routine then 'sub-samples' the generated data, and the undesirable 'piecewise look' in the final plot therefore is still produced.
3. A better approach is to use the semilogplot() command. Maple's adaptive plot routines are now clever enough to produce similar numbers of points with each decade of the independent variable. Obviously using semilogplot() rather than odeplot(), measn that the output of the dsolve() command has to be accessed somewhat differently. However this is no both fast (and accurate!)

odePlot.mw

in Maple 2017.3 on 64-bit Win7, I get

Optimization:-LPSolve(750*xp+250*xq-50000*xifq,maximize=true,
{20*xp+30*xq<=2600,
30*xp+20*xq<=1000,
20*xp+30*xq<=2000,
xq<=99999*xifq,
xq>=20*xifq,
xp<=99999*xifp,
xp>=30*xifp,
xp>=0,
xq>=0},
integervariables={xp,xq},
binaryvariables={xifp,xifq}
);
         [24750, [xifp = 1, xifq = 0, xp = 33, xq = 0]]

However you should be aware that there has been significant recent discussion on problems with LPSolve(), see here

https://www.mapleprimes.com/questions/223522-LPSolve-Doesnt-Register-Constraint

Error message says (my highlighting)

Error, (in pdsolve/numeric/process_PDEs) can only numerically solve PDE with two independent variables, got {t, x, y}

Help page says

PDEsys:   single or set or list of time-dependent partial differential equations in two independent variables

You have three independent variables

When you state that you are trying to solve a " binary problem " do you mean that the variables x__1, x__2,..x__5 can only assume the integer values 0,1? If so defining them with a range 0..1 is completely wrong, because this means that these variable can take on any (floating point) value between 0 and 1

On the other hand the code

restart;
Optimization:-LPSolve( 3*x__1+14*x__2+18*x__3+6*x__4+2*x__5,
                       {3*x__1+5*x__2+6*x__3+2*x__4+x__5 <= 10},
                       binaryvariables=[x__1, x__2, x__3, x__4, x__5],
                       maximize = true
                      );

whihc restricts x__1, x__2, etc to binary values, returns

[26, [x__1 = 0, x__2 = 0, x__3 = 1, x__4 = 1, x__5 = 1]]

in both Maple2016.2 and Maple 2017.3

restart;
eq:= -sin(alpha)*(sin(theta1)*cos(theta)-cos(theta1)*sin(theta));
thaw(combine(subs( sin(alpha)=freeze(sin(alpha)),eq), trig));

The freeze(), thaw() combination is very useful for restricting simplifications/combinations

because now all you have to do is pass the reruired function to be plotted to a procedure, which handles data generation, plot titles, etc, etc

Basic procedure could be modified for different ranges etc etc, just by passing/using additional arguments

odeplot2.mw

using the 'listprocedure' option to the dsolve() command. As in the attached.

odeplot.mw

By the way I think that your original is not plotting what you think you are plotting!

There is one "simple" syntax problem, which was confusing things.

I agree with you that, for the fourth set of constraints, the definition does not  correspond to the subsequently listed equations. I cannot "see" a simple typo in the definition which would produce the listed equations. Maybe there is such a typo - maybe not, just syaing that I can't figure out what it might be.

Using the definition (ie the first group defining 'c4' in your worksheet), the optimiser produces an answer of 4 bins.

Using the second definition of 'c4' in your worksheet (ie matching the corresponding equations in the pdf), the optimser now returns 2 bins (which your original worksheet did not - presumably the impact of the syntax problem mentioned above)

Even this 2-bin solution is not exactly the same as that given in the pdf. 'q', 'x', 'y', values correspond, but 'z' values do not. NB it is not obvious to me that the solution to this problem is necessarily unique, so the solution obtained *may* be correct

Check the attached - my amendments/annotations are coloured red

binPack.mw

By default Maple does floating-point calculations to 10 significant digits, although this can be set to (almost) any value by using the 'Digits' environment variable. The maximum value of Digits which can be set and used is obtained from kernelopts(maxdigits);, which on my machine (64-bit Windows 7, running Maple 2017.3), returns 38654705646 - Now that is a lot of digits!

Most(?) calculators, by default operate with 12 significant digits (by default)  - so 2 more than the Maple default.

However when I try your specific example, for any Digits setting from 4 to 16, I don't see the "rounding error" which you report. See th attached. Possibly your issue is version (or OS?) dependent

Download acc.mw

when

fsolve({A1, A2}, {C1, X});

returns

{C1 = 1.802618458, X = 3.656793458}

although you have to reorganise your input data into a usable format.

Should the first entry for 'B', 'y=10', be [ 4, 10, 2.0] rather than [ 4, 5, 2.0] ?

See the attached for plots

doPlot.mw