If I understand your question, you want to solve the boundary problem for a variety of functions  f(x) . It is advisable to write as a procedure. Only note that this problem may not have a solution. But if it has a solution  u0(x) , it will have the infinite number of solutions  u0(x)+C*sin(x)

restart;
Sol:=F->dsolve({diff(u(x),x,x)+u(x)=F, u(0)=0, u(Pi)=0}):

Example of use:
Sol(-3*sin(2*x));
                                 u(x) = sin(x)*_C1+sin(2*x)

The file below shows your system in the desired form  Sys  (all second derivatives are in the left-hand sides of equations).

Model_Maple2.mw

Your system is complex and has many parameters. I think that such a system can only be solved numerically for specific parameter values and initial conditions. 
Here is my attempt for arbitrary values of parameters and initial conditions:

Model_Maple1.mw

Edit.

Cuts:=proc(n::posint)
local S, R, i, Arcs, P, M, E, PP, d, Eq, Seg, F, Eqs, Dev, a, u, v, r, G, m, j;
global C, IC, L;
uses plottools, plots;
C:=circle(color=blue, thickness=3);
IC:=disk(color="LightYellow");
S[1]:=[-1.2, 1.2]; # The first two points
Seg[1]:=[cos,sin]~(S[1]);  # The first line segment specified by the coordinates of its ends
L[1]:=line(Seg[1][], color=red, thickness=1); 
R:=[]; Eqs:=[];
for i from 1 to n-1 do
R:=sort([op(R),op(S[i])]); # The sorted list of points on the circle
Arcs:=[seq([R[i],R[i+1]], i=1..nops(R)-1), [R[-1],R[1]+2*Pi]]; # The arcs on the circle, formed by the points R
M:=map(t->abs(t[1]-t[2]),Arcs); d:=max(M); m:=nops(Arcs);
select(t->abs(t[1]-t[2])=d, Arcs)[1];
E:=(%[1]+%[2])/2; # The middle of the maximum arc
PP:=map(t->(t[1]+t[2])/2, Arcs);  # The list of points of all the arcs
F:=(x,y,Q)->(x-Q[1,1])*(Q[2,2]-Q[1,2])-(Q[2,1]-Q[1,1])*(y-Q[1,2]);
Eqs:=[op(Eqs), F(x,y,Seg[i])];
Dev:=map(t->[seq(eval(Eqs[i],{x=cos(t),y=sin(t)}),i=1..nops(Eqs))], PP);
a:=[seq(eval(Eqs[i],{x=cos(E),y=sin(E)}),i=1..nops(Eqs))];
for j from 1 to m do
[seq(a[k]*Dev[j,k], k=1..i)];
if `and`(seq(%[k]<0, k=1..i)) then break fi;
od;
u:=0.9*Arcs[j,1]+0.1*Arcs[j,2];  v:=0.1*Arcs[j,1]+0.9*Arcs[j,2];
r:=`if`(u<0 and v>0,rand(0..v-u),rand(u..v));
G:=`if`(u<0 and v>0,r()+u,r());
S[i+1]:=[E,G];
Seg[i+1]:=[cos,sin]~(S[i+1]); 
L[i+1]:=line(Seg[i+1][], color=red, thickness=1);
od;
display(C, IC, seq(L[i],i=1..n), axes=none);
end proc:

Examples of use.

The plot:

Cuts(5);

The animation:

F:=n->plots:-display(C, IC, seq(L[i],i=1..round(n)), axes=none, scaling=constrained):
plots:-animate(plots:-display,['F'('n')], 'n'=0..5, frames=60);

Cuts.mw

Edit. The code of the procedure was edited.

In your code, the parameter name  r  matches the name of the function  r(t) . I changed the name of the parameter and for the plotting I have reduced the range for variable t :

restart;
beta:=0.00013: R:=0.03:
sir:={diff(s(t),t)=-beta*s(t)*i(t), diff(i(t),t)=beta*s(t)*i(t)-R*i(t), diff(r(t),t)=R*i(t)};
DEtools:-DEplot(sir, [s(t),i(t),r(t)], t=0..30, scene=[t,i(t)], [[s(0)=10000, i(0)=5, r(0)=0]]);

in addition to the rlopez's one:

restart;
P := x + y + a*x^2 + b*x^3 + c*x*y + d*x^2*y + e*x*y^2 + y^2:
select(i->degree(i, [x,y])=2, P);     # or
map(i->`if`(degree(i, [x,y])=2, i, 0), P);

                        a*x^2+c*x*y+y^2
                        a*x^2+c*x*y+y^2
 

U := proc (p, q, R)
local a;
Digits:=5;
a := arccos(1/2*(p^2+R^2-q^2)/(p*R));
alpharad = a, alphadeg = 180*a/Pi; 
end proc:

Example:

select(t->rhs(t)=0,  [x[`001`]=0, x[101]=1, x[201]=0, x[301]=1]);

                                     [x[`001`] = 0, x[201] = 0]

 or

remove(t->rhs(t)=0,  [x[`001`]=0, x[101]=1, x[201]=0, x[301]=1]);

                                        [x[101]=1, x[301]=1]

Instead of   i, j -> ...   should be   (i, j) -> ...   (parentheses are missed)

See corrected file:

h_2.5_eq4_new.mw

XY:=[[1,2.58], [3,8.57], [5,15.00], [10,30.10], [15,45.00]]:
CurveFitting:-LeastSquares(XY, x);
plot([XY, %], x=0..16, style=[point,line], symbolsize=15);

1. For exact solution of a linear system use  LinearAlgebra:-LinearSolve  or  solve   command.  

2. If the system has more equations than unknowns, it is usually inconsistent, but it can be solved approximately using   LinearAlgebra:-LeastSquares  command .

3. See help on these commands.

Actually we have geometric series with  q=sin(k*x/2) :

eval(sum(q^(2*n-2), n=0..N), q=sin(k*x/2));  # Finite sum
S:=sum(q^(2*n-2), n=0..infinity) assuming abs(q)<1;  # Infinite sum
eval(S, q=sin(k*x/2));

Edit. If the coefficients of the sum are not equal to  1 , the change  q=sin(k*x/2)   may also be useful: first, by simplifying the expression of  q , and then doing the reverse change.

content(4-2*a);
primpart(4-2*a);
%%*``(%);

A trick is needed for the second example:

P:=2*x^2 - 2*x + sqrt(2);
sqrt(2)*``(expand(P/sqrt(2)));
map(`^`,  %, 2);

Edit. The first example can be written shorter:

content(4-2*a, a, p)*``(p);

                                               2*(2-a)

Here is a solution to your problem on a simple example:

A:=Matrix([[a1*x+b1*y+c1*z],[a2*x+b2*y+c2*z],[a3*x+b3*y+c3*z]]);
LinearAlgebra:-GenerateMatrix(convert(A,list), [x,y,z])[1];  # The desired matrix
%.<x,y,z>;  # Verification

deqv := m*v(s)*(diff(v(s), s)) = m*g-k*v(s)^2;
solv := dsolve({deqv, v(0) = v0}, v(s))  assuming m>0, k>0, g>0, v0>0;