In this system, each equation contains only one unknown. Therefore each unknown can be found by the same formula with two parameters  a=sum(a[i], i=1..n)  and  b :

solve(a*x+x^3=b, x);

See  Wiki   for details.

Below we find the general and a particular solution of the equation for specific parameters:

restart;
Eq:=P=i(t)^2*r+diff(1/2*l(c)*i(t)^2, t)+1/2*i(t)^2*diff(l(c), c)*w:
dsolve(Eq, i(t));
Eq1:=eval(Eq, {P=2, r=3, l(c)=sin(c), w=1});
dsolve({eval(Eq1, c=Pi/3), i(0)=5}, i(t));

The recursive procedure  Prefix  automates the problem for 3 operations  `*`, `+`, `^` :

Prefix:=proc(Expr)
if type(Expr, `*`) then return func1(Prefix~([op(Expr)])[ ]) else
if type(Expr, `+`) then return func2(Prefix~([op(Expr)])[ ]) else
if type(Expr, `^`) then return func3(Prefix~([op(Expr)])[ ]) else
Expr  fi; fi; fi;
end proc:

Examples of use.

The original example:

Prefix(a^2*b+c);
                                                   

A more complicated example:

Expr:=(a*b*c+d^2)*(a+b*c^(-1))^3;
A:=Prefix(Expr);
subs({func1=`*`,func2=`+`,func3=`^`}, A);
expand(%);  # Check
expand(Expr);

Addition. The procedure also works for division because

op(a/b);
whattype(a/b);

                                          a, b^(-1)
                                                *

You mean the prefix notation? In Maple it is implemented as (for your example)

restart;
L:=[a, b, c]:
`+`(`*`(L[1]$2, L[2]), L[3]);

                                                       a^2 * b + c
                               

Addition. Prefix notation is useful if you want to apply an operation to a large or indeterminate (in advance) number of operands.

A simple example - find the sum of the first 100 natural numbers:

`+`($ 1..100);
                                           5050

`if`(modp((n-1)!=-1, n^2), n, NULL) $ n=1..1000;

                                                                      1, 5, 13, 563

seq(`if`(modp((n-1)!=-1, n^2), n, NULL), n=1..1000);   # by  seq  command

                                                                      1, 5, 13, 563

for n from 1 to 1000 do   # by a  for  loop
if modp((n-1)!=-1, n^2) then print(n) fi;
od;

                                                                            1
                                                                            5
                                                                           13
                                                                          563
 

Edit.

3_ways.mw

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