Your sequence can be defined by recurrent formula or  procedure-function.

By recurrent formula:

By procedure:

a:=proc(n)

local b, i;

b[1]:=1; b[2]:=2;

for i from 3 to n do

if type(i, odd) then b[i]:=b[i-2]*b[i-1] else b[i]:=b[i-3]+b[i-2] fi;

od;

b[n];

end proc: 

Example:

seq(a(n), n=1..20);

1, 2, 2, 3, 6, 5, 30, 11, 330, 41, 13530, 371, 5019630, 13901, 69777876630, 5033531, 351229105131280530, 69782910161, 24509789089304573335878465330, 351229174914190691

restart;

with(plots):

A:=contourplot(x^3-y^2,x=-3..3,y=-3..2, thickness=2, axes=boxed, contours=[seq(-4..5)]):

B:=textplot([seq([surd(0.7^2+i,3),0.7, i], i=-4..5)], font=[TIMES,ITALIC,12], color=brown, align={left, above}):

display(A,B);

There is no  labelledcontourplot  command in Maple.

plots[contourplot](x^3-y^2,x=-3..3,y=-3..2, thickness=2, axes=boxed, contours=10);

3 options are working:

L:=[sum, value, plot]:

for i in L do

Student[Calculus1][ApproximateInt](sin(Pi*x), x=0..1, functionoptions=[filled=true, color=pink], output=i);

od;

If "to reduce"  means reduction in the number of lines of code, then possible variant

simplify(solve(map(t->t<2, [solve(x^4 -(3*m+2)*x^2 + 3*m+3=0, x)]), m));

                              {2/3*2^(1/2) <= m, m < 11/9}

To improve the quality of plot use  numpoints  option, for example

DEplot({de1, de2}, [x(t), y(t)], t = 0 .. 5, [[x(0) = 10, y(0) = 0]], stepsize = 0.5,

color = black, linecolor = blue, numpoints=5000);

For the second and third plots arrows are not drawn because the equations are non-autonomous.

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);