@ReactionUra Usually, when fsolve returns unevaluated, it means that there's no solution within the specified ranges. However, sometimes it means that it couldn't get the solution to sufficient precision. That's the case here. By reducing the requested precision to 7 digits, I got an all-nonnegative solution:

Digits:= 7:
fsolve(
    {eq||(1..12)},
    {ionic, n_Cl, n_H, n_HCl, n_Na, n_NaOH, n_OH, m, u_Cl, u_H, u_Na, u_OH}=~ 0..infinity    
);
  /                                                     
 { ionic = 48.98702, m = 0.003277555, n_Cl = 0.1750121, 
  \                                                     

   n_H = 0.05026130, n_HCl = 0.00001712527, n_Na = 0.09573009, 

                                  -13                    
   n_NaOH = 0., n_OH = 1.128623 10   , u_Cl = 0.6624931, 

                                                   \ 
   u_H = 963.1730, u_Na = 137.5875, u_OH = 459.2063 }
                                                   / 

@janhardo Okay, I added some exposition to the worksheet. Let me know if you can understand it now. Unfortunately, MaplePrimes won't display it, but you can download it. By the way, the derivatives that I used are not implicit.

Maple Worksheet - Error

Failed to load the worksheet /maplenet/convert/PointPlaneDistance.mw .
 

Download PointPlaneDistance.mw

@ReactionUra 

Use fsolve instead of solve:

fsolve({eq||(1..12)});
res := {ionic = -2.247127086, m = 0., n_Cl = -3.680938534, n_H = -0.2101029014, n_HCl = -7.502816177*10^7, n_Na = 4.769998672, n_NaOH = -4.769878672, n_OH = 8.240834305, u_Cl = -4.874160233, u_H = 3.688110962, u_Na = -2.974047722, u_OH = -4.110097263}

I suspect that some of these values are outside of allowed ranges. In particular, I suspect that the n_ variables should be nonnegative. You can specify the allowed ranges like this:

fsolve(
    {eq||(1..12)},
    ({n_Cl, n_H, n_HCl, n_Na, n_NaOH, n_OH}=~ 0..infinity)
    union
    ({ionic, m, u_Cl, u_H, u_Na, u_OH}=~ -infinity..infinity)
);
{ionic = -2.225127530, m = 0., n_Cl = 3.159530733, n_H = 5.462749423, n_HCl = 2.054664202*10^7, n_Na = 7.562418530, n_NaOH = 3.102492997*10^6, n_OH = 9.037400983, u_Cl = -4.874160233, u_H = 3.641064071, u_Na = 1.655806581, u_OH = -4.110097263}

@vv In this particular case, the results are practically the same without using parametric or real. I mention this because I think that the OP is using a very old version of Maple where the parametric option may work differently, if it works at all.

@ReactionUra Which are the 13 variables? Using Maple 2020, I solved it almost immediately, but I let it pick the variables to solve for. I can't investigate further without knowing the variables.

Also note that there are two instances of exp^(...) that need to be corrected.

I can't tell what's going on without seeing how you've used Units in the preceding code. But this should give you the correct result:

max~(0, ResLoadcombination_unitfree)

Does that work for you?

Please Post this as a Reply to your earlier Question.

If you're trying to find a value of L such that y(2) = 0, it's done already. That value is L = 2.717402017. There's no need to use bisection as a root-finding method; fsolve takes care of that.

If you're going to download the OP's worksheet, make sure to turn off autosave first. If you don't, the first autosave will never finish, and you'll be forced to kill your entire Maple session. To turn off autosave, go to Tools => Options => General, uncheck the last item, and Apply to Session.

@janhardo Okay, I see my mistake of a lack of generality. Your solution can be obtained from pdsolve by

pdsolve(diff(u(x,y), x, y) = -diff(u(x,y), x));
             u(x, y) = _F1(y) + exp(-y)*_F2(x)

@janhardo What long calculation? I get the solution _C1*x*exp(-y) + _C2.

@acer Acer, would you please explain the purpose of

a:= 'a':  b:= 'b':  c:= 'c':

that appears before the gc() that appears before timed code? My guess is that its purpose is to force garbage collection of %, %%, and %%%.

Also, I'm guessing from the way that you recoded ZigZag that compiled procedures can't use op(procname). I used that to avoid the extra procedure layer (i,j)-> ZigZag(i,j,n).

@Rouben Rostamian  Thanks. When you edit a Reply, there's a pull-down at the bottom that says Parent. That's how you move Replies to different Answers. 

@Rouben Rostamian  Here is the inverse of the above procedure, i.e., it returns the indices given the entry:

ZigZagInv:= (k::posint)-> 
local n:= op(procname);
    if k > n*(n+1)/2 then n+1 -~ thisproc(n^2+1-k)
    else
        local d:= ceil((sqrt(1+8*k)-1)/2), r:= k - (d-1)*d/2;
        [r, d+1-r][(-1)^d*[1,2]]
    fi
:
#Example usage
M:= Matrix(8$2):
for k to 64 do M[ZigZagInv[8](k)[]]:= k od:
M;

@acer The design that units are names rather than strings (or records or objects) seems to me to be inherently flawed.

@Sayed From my code above, the set of coefficients alone is simply {C}.