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]

An example:

A:=plot(sin(x), color=red,legend="Sin"):

B:=plot(cos(x), color=blue, legend="Cos"):

C:=plot(tan(x), color=green, legend="Tan", discont=true):

E:=plot([seq([Pi/2+Pi*k, t, t=-3..3], k=-1..2)],color=black, linestyle=2, legend=["Asymptote"$4]):

plots[display](A,B,C,E, thickness=2,view=[-Pi/2..5*Pi/2, -3..3], scaling=constrained);

Corrections in the last line and result:

plots[implicitplot]('IC(2, x, y, 1, 1)', x = .1 .. .9, y = -1 .. 1, filled = true, numpoints = 10);

#  or  plots[implicitplot]((x,y)-> IC(2, x, y, 1, 1), .1 .. .9, -1 .. 1, filled = true, numpoints = 10);