Since solutions of  your system  [Er1, Ez1]  are symmetric respect to the line  rho1=0 , then we assume  rho1>0 . It is easy to check that the second equation has a solution   {z1=R/2, rho1 - any number} ,  that will be unique for  R>0  and approximately  R<=1.2 .For such values ​​of  R  the second unknown  rho1  can be obtained from the solution of the first equation for  z1=R/2 . This solution will be called trivial. The blue line on the Carl Love's plot  in the range  0<R<=1.2  is exactly a straight line with equation z1=R/2 .

If  R> 1.2, then except the trivial solution there are two non-trivial solutions.  Carl Love's code finds smaller solution  z1[1]<R/2 . Since  non-trivial solutions are symmetric respect to the line  z1=R/2 , then the second non-trivial solution be  z1[2]=R - z1[1] . 

All of the above conclusions are clearly seen from the plots of the equations of the system. Plots are made for specific  R  in the range  0.2<=R<=2  with step=0.2 . 

Branching_of_solutio.mw

I just copied your code into a new worksheet. Everything works!

TriDiagMatr.mw

linsolve is the command in  linalg  package. First you need to download the  linalg  package, and then use  linsolve  command. 

Example:

with(linalg):

A := matrix( [[1,2],[1,3]] ):

b := vector( [1,-2]):

linsolve(A, b);

         [7, -3]

Two examples.

1) For  m=1:

A:=TriDiagMatr(1):

x:=Vector([seq(a[i], i=-1..2)]):

b:=Vector([1,7,5,10]):

solve({seq((A.x)[i]=b[i], i=1..4)});

2) If  m> 1 , then the symbolic solution is too cumbersome, so rather than  solve  it is better to use  fsolve

A:=TriDiagMatr(5):

x:=Vector([seq(a[i], i=-1..6)]):

b:=Vector([1,7,9,11,3,5,10,13]):

fsolve({seq((A.x)[i]=b[i], i=1..8)});

The procedure  TriDiagMatr  allows to build your three-diagonal matrix for any  m>=1.

TriDiagMatr:= proc(m)

local L, i, B, s;

L:=Array(0..m, [seq(((a[i-1]+4*a[i]+a[i+1])/6)^2, i=0..m)]);

for i from 0 to m do

B[i-1,i]:=6*i+24+L[i];

B[i+1,i]:=6*i-24+L[i];

B[i,i]:=6*i+L[i];

od;

s := {seq(seq((i,j) = B[j-2,i-2], j=i-1..i+1 ), i = 2 .. m+2), (1, 1) = 1, (1, 2) = 4, (1, 3) = 1, (m+3, m+1) = 1, (m+3, m+2) = 4, (m+3, m+3) = 1};

Matrix(m+3, s);

end proc:

Example for m=3:

TriDiagMatr(3);

You can see the process of solving step by step as follows:

with(Student[Calculus1]):

A:=Int(1/(x^2+2), x):

ShowSolution(A);

After the asymptotes found you build them together with the  original function on the same plot.

Example.

with(Student[Calculus1]):

Asymptotes(1/(x-3)+2*x, x);

plot([1/(x-3)+2*x, 2*x, [3, t, t = -10 .. 20]], x = -2 .. 7, -10 .. 20, color = [red, black, black], thickness = [2, 1, 1], linestyle = [1, 2, 2], discont = true);

P:=proc(`%`)

end proc;

Error, cannot use environment variable `%` as a parameter in procedure P

In one command line:

plot([2, 5, seq(k, k = 2.05 .. 4.95, 0.1)], x = -2 .. 6, -2 .. 6, color = [red, "DarkGreen", "LightBlue" $ 30], thickness = [6, 5 $ 31], gridlines = true);

In your example incorrect syntax. The name  f(x)   should be assigned  to  the expression  x^2+1 .  Then all works. The maximum value of this expression can be found exactly (not approximately) by  maximize  command:

f(x) := x^2+1;

Xmax:=maximize(f(x), x=0..10);

plot(f(x)/Xmax, x=0..10);

For  b=0.1  the problem can easily be solved exactly. It is interesting to compare exact and numerical solutions (with numeric option).

b:=0.1:  L:=[seq(x, x=0...1, 1/10)];

dsolve({(x+b*y(x))*diff(y(x), x) + y(x) =0, y(1)=1});

seq(eval(rhs(%),x=L[i]), i=1..nops(L));

evalf[18]([%]);

Let diagonally for  i=3 .. m+3  as in your example are the same numbers (elements of the list  L ).

Matr1 := proc(m, L) 

local s;

s := {seq(seq((i, j) = L[j-i+3], j = i-2..i+2), i = 3 .. m+3), (1, 2) = 1, (1, 3) = 4, (1, 4) = 1, (2, 2) = -1, (2, 4) = 1, (m+4, m+2) = -1, (m+4, m+4) = 1, (m+5, m+2) = 1, (m+5, m+3) = 4, (m+5, m+4) = 1};

Matrix(m+5, s);

end proc:

Example for  m=5:

Matr1(5, [6, 2*x, 3, 4, 9]);

According to the problem, the roots must be different. The result may be easily verified graphically. And the initial problem itself also can be solved graphically.

sol:=[m, -3*m+2, -m-3, 2*m+1]:

solve({max(sol) > 0, mul(j, j = sol) > 0, sol[1] <> sol[2], sol[1] <> sol[4], sol[3] <> sol[2], sol[3] <> sol[4], min(sol) < 0}, m);

plot([m, -3*m+2, -m-3, 2*m+1], m = -5 .. 5, -8 .. 10, color = [red, blue, green, yellow], thickness=2);

You seem to have missed the element  a[4,2] .

The procedure  Matr  solves your problem, for any  m>=0 .

Matr := proc(m)

local s;

s := {seq(seq((i, j) = a[j-3, i-3], j = i-2 .. i+2), i = 3 .. m+3), (1, 2) = 1, (1, 3) = 4, (1, 4) = 1, (2, 2) = -1, (2, 4) = 1, (m+4, m+2) = -1, (m+4, m+4) = 1, (m+5, m+2) = 1, (m+5, m+3) = 4, (m+5, m+4) = 1};

Matrix(m+5, s);

end proc:

Your example:

Matr(3);

Another example:

interface(rtablesize=20):
Matr(7);

Use the command  Student[VectorCalculus][LineInt]