From the link   webpage  you specified, you missed the beginning of the code, in which all parameters  (a, b, h, n, f)  were set. Here is the full code for your version (you have changed only the parameter  n ):

Download Adjust.mw

The equation  DEF_Profit1  in addition to the unknown  S  contains 3 parameters  beta, s, C. If you specify the values of these parameters, you can easily find all the roots, for example:

sol1 := solve(eval(DEF_Profit1, [beta = 3, s = 0.1, C = 2]), S);
                        sol1 := 0.004990229958, 29.61213284, -3.574031930, -26.04309114

By default, axes=normal  option works. So just remove  axes=box  option from your code:

plot(sin(x), x=0..2*Pi);

Addition.
 If you want the boxed effect and that the tickmarks are on the line  y=0 , then do so:

plots:-display(plot(sin(x), x=0..2*Pi), plot([[0,-1],[2*Pi,-1],[2*Pi,1],[0,1]], color=black));
                

Edit.         

If you do not have additional information about the parameters  chi, lambda1, lambda2, phi, psi, rho, you can simply use  simplify  command:

with(LinearAlgebra):

A:=Matrix([[phi,conjugate(psi),chi,conjugate(rho)],[psi,-conjugate(phi),rho,-conjugate(chi)],[lambda1*phi,conjugate(lambda1)*conjugate(psi),lambda2*chi,conjugate(lambda2)*conjugate(rho)],[lambda1*psi,-conjugate(lambda1)*conjugate(phi),lambda2*rho,-conjugate(lambda2)*conjugate(chi)]]);
indets(A, name);
Determinant(A);
simplify(%);
                           

Try  colorscheme  option for curves. An example:

restart;
plots:-spacecurve([cos(t), sin(t), 0.3*t], t=0..4*Pi, colorscheme=[blue, green], thickness=3, axes=normal);
                            

Unfortunately, this does not work when solving differential equations with numeric option (plotting by  plots:-odeplot  command). In this case, use a workaround as below (the same curve as above):

restart;
sol:=dsolve({diff(x(t),t)=-sin(t), diff(y(t),t)=cos(t),diff(z(t),t)=0.3, x(0)=1,y(0)=0,z(0)=0}, numeric):
plots:-pointplot3d([seq(eval([x(t),y(t),z(t)],sol(i)), i = 0..12.56,0.2)], colorscheme = [blue, green], style = line, thickness=3, axes=normal);
 

You can specify  g(x)  as a piecewise-function and then use in any calculations, for example:

f(x):=x+1;
g(x):=piecewise(x<>2,(x+1)*(x-2)/(x-2), x=2,undefined);
simplify(f(x)*(g(x)+1));
limit(g(x), x=2);
diff(g(x), x);
int(g(x), x=1..3);
discont(g(x), x);
                                   

With  for  loop:

restart;
r:=rand(1..100):
k:=0:
for i from 1 to 20 do
L[i]:=r();
if type(L[i], even) then k:=k+1; L1[k]:=L[i] fi;
od:
L:=convert(L, list);
L1:=convert(L1, list);
`+`(L1[])/nops(L1);

It is easier to do this in another way:

restart;
r:=rand(1..100):
L:=[seq(r(), i=1..20)];
L1:=select(t->type(t, even), L);
`+`(L1[])/nops(L1);

If you want new random numbers to be generated with each call, then insert the  line  randomize():  after restart;

Programmatically, I found only a method using  plots:-textplot  command:

plots:-textplot([[-0.1,0, max,font=[times,bold,20]], [0.1,0, f(x),font=[times,20]], [0,-0.1, x in [2,3]]], size=[140,60], axes=none);
                                                              

1. Conversion to  Matrix  can be easier (see the code below).

2. I agree that this piece "A particular solution of the system can be then written in the form   F(t) = S(t) C0 + P(t)  where  C0 is  n by 1 and  F(0) = C0 + P(0)." of Help is incorrect.

3. The constants  C1  and  C2  for the particular solution with the specified initial conditions  x0  can be obtained as follows (first we find the general solution  x ):

restart;
A:=Matrix(2, 2, [1, 4, 5, 1]);
x0:=Vector([2,1]); #initial condition
sol := DETools:-matrixDE(A, t);
M:=convert(sol[1], Matrix);  #convert to Matrix
x:=M.<C1,C2>;  # C1 and C2 are arbitrary constants
solve(eval(convert(x=~x0,list),t=0), [C1,C2])[];
simplify(eval(x, [%[],t=0])); # The check

dice := proc()
    local r;
    r:=rand(1..6);
    'r'()$2;
end proc:

Example of use:
seq([dice()], i=1..10);
             [6, 5], [3, 1], [5, 2], [3, 2], [2, 4], [3, 3], [1, 2], [5, 4], [5, 6], [2, 3]

The code below finds 2 chains:

chosen1 := [[1,2],[1,20],[3,4]]:
chosen2 := [[2,3],[20,3],[3,4]]:
chosen3 := [[3,4],[5,7]]:
chosen4 := [[4,5],[5,6]]:
chosen5 := [[5,6],[7,9]]:
seq(seq(seq(seq(seq(`if`(i[2]=j[1] and j[2]=k[1] and k[2]=l[1] and l[2]=m[1], [i,j,k,l,m], NULL), m=chosen5), l=chosen4), k=chosen3), j=chosen2), i= chosen1);
                       [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6]],  [[1, 20], [20, 3], [3, 4], [4, 5], [5, 6]]
 

Addition.  If there are many source lists, then the process of building nested sequences can be automated using the procedure  NestedSeq  that can be found in this thread.

because

evalc(cos(I*x));  # This assumes that x is real
convert(%, exp);
                                                cosh(x)
                                     (1/2)*exp(x)+(1/2)*exp(-x)

Use  extrema  command for this (gives min and max). Optimization:-Minimize  command gives the same result.

An example:

Download min.mw

See help on  CurveFitting:-PolynomialInterpolation  command.

A:=[[1,2],[5,7]]:
B:=[[1,2]]:
(convert(A,set) intersect convert(B,set))[];

# or

A:=[[1,2],[5,7]]:
B:=[[1,2]]:
({A[]} intersect {B[]})[];