I'm always skeptical of assuming some_property without specifying variables. If you change it to assuming x::real, then it works. Otherwise, some problem occurs when it tries to also assume that right is real. I don't know why that only happens when you use quotes.

I assume that you know that the assumption real is redundant when using left or right, but for whatever reason you want to use it anyway.

Another solution is to make right protected. (I guess that assumptions on unspecified variables don't apply to protected variables.)

protect(right);
...rest of your code...

The mathematical constant is spelled Pi (capital P); pi (lowercase) is just another variable.

To convert the units, do 

convert(Rac, units, ohm);

The displayed form will use the capital Omega that you expect.

Here's a procedure for it. Since I suspect that you only want to exclude the integral if the special function depends on the variable of integration, I included an optional second parameter for a variable name or set of variable names. The procedure is constructed with subs so that the calls to FunctionAdvisor and the construction of the set are only done once.

has_special_math_function:= subs(
    _F= {(op@FunctionAdvisor)~(FunctionAdvisor("class_members", "quiet"), "quiet")[]}
        minus {FunctionAdvisor("elementary", "quiet")[]},
    proc(e::algebraic, x::{name, set(name)}:= {})
        hastype(e, And(specfunc(_F), dependent(x)))
    end proc
);
has_special_math_function := proc (e::algebraic, x::{name, set(name)} := {})
    hastype(e, And(specfunc({AiryAi, AiryBi, 
        AngerJ, BellB, BesselI, BesselJ, BesselK, BesselY, Beta, 
        Chi, Ci, Dirac, Ei, GAMMA, HeunB, HeunC, HeunD, HeunG, 
        HeunT, Im, JacobiP, KummerM, KummerU, Li, MeijerG, Psi, Re, 
        Shi, Si, Ssi, StruveH, StruveL, WeberE, Zeta, abs, dawson, 
        dilog, erf, erfc, erfi, euler, fourier, hankel, hilbert, 
        laplace, lnGAMMA, max, mellin, min, polylog, signum, 
        unwindK, AppellF1, AppellF2, AppellF3, AppellF4, ChebyshevT, 
        ChebyshevU, CompleteBellB, CoulombF, CylinderD, CylinderU, 
        CylinderV, EllipticCE, EllipticCK, EllipticCPi, EllipticE, 
        EllipticF, EllipticK, EllipticModulus, EllipticNome, 
        EllipticPi, FresnelC, FresnelS, Fresnelf, Fresnelg, 
        GaussAGM, GegenbauerC, HankelH1, HankelH2, Heaviside, 
        HermiteH, HeunBPrime, HeunCPrime, HeunDPrime, HeunGPrime, 
        HeunTPrime, IncompleteBellB, InverseJacobiAM, InverseJacobiCD,
        InverseJacobiCN, InverseJacobiCS, InverseJacobiDC, 
        InverseJacobiDN, InverseJacobiDS, InverseJacobiNC, 
        InverseJacobiND, InverseJacobiNS, InverseJacobiSC, 
        InverseJacobiSD, InverseJacobiSN, JacobiAM, JacobiCD, 
        JacobiCN, JacobiCS, JacobiDC, JacobiDN, JacobiDS, JacobiNC, 
        JacobiND, JacobiNS, JacobiSC, JacobiSD, JacobiSN, 
        JacobiTheta1, JacobiTheta2, JacobiTheta3, JacobiTheta4, 
        JacobiZeta, KelvinBei, KelvinBer, KelvinHei, KelvinHer, 
        KelvinKei, KelvinKer, LaguerreL, LambertW, LegendreP, 
        LegendreQ, LerchPhi, LommelS1, LommelS2, MathieuA, MathieuB, 
        MathieuC, MathieuCE, MathieuCEPrime, MathieuCPrime, 
        MathieuExponent, MathieuFloquet, MathieuS, MathieuSE, 
        MathieuSEPrime, MathieuSPrime, MultiPolylog, MultiZeta, 
        NielsenPolylog, SphericalY, Stirling1, Stirling2, 
        WeierstrassP, WeierstrassZeta, WhittakerM, WhittakerW, 
        Wrightomega, argument, bernoulli, binomial, conjugate, 
        doublefactorial, factorial, fouriercos, fouriersin, 
        harmonic, hypergeom, invfourier, invhilbert, invlaplace, 
        invmellin, multinomial, piecewise, pochhammer, GeneralizedPolylog, 
        MathieuFloquetPrime, WeierstrassPPrime, WeierstrassSigma}), 
   dependent(x))) 
end proc

You may want to exclude some functions such as signum, abs, Re, Im, argument, etc. from that set.

Setting the value of _Envsignum0 only affects signum(0), not signum(a) for arbitrary a.

FWIW, you could do 

limit(diff(y*fy, y), y= 0, right) assuming a > 0

to make signum(a) = 1. For that matter, you could do

eval(limit(diff(y*fy, y), y= 0, right), [signum= 1, csgn= 1])

although it's not clear to me why you don't want signum in the result.

The command 

series(x/(tan(x) - sin(x)), x= 0)

shows that f(x) is approximately equal to 2/x^2 for x close to 0. This suggests that g(x) = 1/x^2 might be a reasonable choice.

In procedure SOR, change x_0:= x to

x_0[]:= x

This overwrites the vector's elements rather than its pointer.

Like this:

plots:-animatte(plot, [[[x1, z1], [x2, z2](t), [x3, z3](t)]], t= 0..T, frames= n);

where T is positive and n a positive integer.

If you want to create a dynamically sized array row-by-row (rather than element-by-element) and then convert it to a matrix, it can be done like this:

A:= Array(1..0):
for i to 9 do A,= (i, i+1, i+2) od:
B:= ArrayTools:-Alias(A, [9,3], C_order);

This is an efficient conversion requiring no intermediate lists or other structures. It doesn't even need to copy the data; it just creates a new pointer to the original array with a different indexing scheme.

The output of plottools:-transform, which you've assigned to t1, is a procedure, so you need to use procedure-calling syntax to use it. Change t1{p1, p2, p3} to

t1~([p1, p2, p3])

Or, you can use t1 to operate on the combined plot like this:

t1(plots:-display(p1, p2, p3))

If you're performing a chain of substitutions, it doesn't make sense to me to introduce intermediate variables that are both

1. not wanted in the final result,
2. not eliminated by some part of the chain.

That's the case with your f0, and it's what causes your troubles. Just remove f0 from equ5.

You've also made a logical error in equ4: `w__k,norm` should be w__k because `w__k,norm` is already eliminated by equ1.

The command subs is inappropriate for doing sophisticated algebra; a better choice is simplify with side relations (see ?simplify,siderels).

After making those corrections, this will do the conversion that you want:

expand(
    simplify(
        errfunc, combinedSet,
        #Since combinedSet has 7 equations, it's best to choose 7 variables to eliminate, but
        #there's some flexibility allowed in which you choose and in which order you put them. 
        [Z__k, Z0, R__x, C__x, L__x, w__k, z__k]
    )
);

Let's say that your original sequence of solutions is named sols. Then do

solsnn:= select(type, {sols}, set(name=nonnegative))[];

The Answer by @mmcdara is very educational. You can learn a lot about systems of equations and about solving them with Maple from it, and you should study it. Nonetheless, there is an easier way which isn't so educational.

restart:
(ia, va, ib, vb):= Vector~(4, symbol=~ [Ia, Va, Ib, Vb])[]: 
(ya, yb):= Matrix~((4,4), symbol=~ [Ya, Yb])[]:
ConnectCondition:= 
    (Ia[3], Ia[4], Va[3], Va[4])=~ (-Ib[1], -Ib[2], Vb[1], Vb[2]) 
:
Y11Condition:= (Va[2], Vb[3], Vb[4])=~ 0:
Eqs:= eval(
    {seq(ya.va - ia), seq(yb.vb - ib)}, 
    {ConnectCondition, Y11Condition}
):
Svars:= seq~({ia, ib, va, vb}) intersect indets(Eqs) minus {Va[1]}: 
S:= subsindets(
    simplify(solve(Eqs, Svars)),
    patindex(anything),
    n-> cat(op(0,n), op(n))
):
lprint(Ia1 = eval(Ia1, S));
Ia1 = (((Ya44+Yb22)*Ya33+(-Yb21-Ya43)*Ya34+Ya44*Yb11+Yb11*Yb22-Yb12*Yb21-Ya43*
Yb12)*Ya11-Ya14*Ya33*Ya41+Ya13*Ya34*Ya41+Ya14*Ya31*Ya43-Ya13*Ya31*Ya44-Ya14*
Ya41*Yb11+Ya13*Ya41*Yb12+Ya31*(-Ya13*Yb22+Ya14*Yb21))*Va1/((Ya44+Yb22)*Ya33+(-
Yb21-Ya43)*Ya34+Ya44*Yb11+Yb11*Yb22-Yb12*Yb21-Ya43*Yb12)

The & is the beginning of an HTML code for a special character; `>` is the "greater than" symbol >; `"` is the double-quotation symbol ". These are not part of ordinary Maple code (unless entered as a part of a quoted string). You must be copying from an HTML source. In actual Maple code, those 4 lines are

for m from n by -2 while n > 2 do
if errr > perrr or errrs > perrrs then
if derrr > 10^10 then
printf(" \134n");

I don't know what \134n in the last command means. Perhaps it's supposed to be just \n, which is the newline character.

All of that (isotherms, 3d plots of Nusselt number, skin friction, etc.) is covered in this Post: https://www.mapleprimes.com/posts/209715-Numerically-Solving-BVPs-That-Have-Many

Like this:

PlotEllipseFromMatrix:= proc(H::Matrix(2, realcons), C::Vector(2, realcons))    
local x, y, X:= <x,y>-C;
    if LinearAlgebra:-IsDefinite(H) = true then 
        algcurves:-plot_real_curve(X^+.H.X - 1, x, y)
    else
        error "not an ellipse, or it can't be determined"
    fi
end proc
:

 PlotEllipseFromMatrix(<3,1;1,1>, <1,2>);