There is another simple way to prevent the immediate calculation of a trigonometric function. If the argument contains  Pi , then just replace  Pi  with  pi . On the output  pi  will look exactly the same as  Pi . When using the  % symbol with the name of a constant or function, they will be displayed in gray, which is not always desirable:

A:=%cos(2*Pi/3)+I*%sin(2*Pi/3);
B:=cos(2*pi/3)+I*sin(2*pi/3);
eval(B, pi=Pi);
                              

Do not use the name and the indexed name with the same symbol in the same document. In your code these are  c  and  c[1], c[2], c[3] . Take a look at what might happen in this case:

restart;
c:=6:
c[1];
                              

I replaced  с  with  С , but then a new error appears. Probably the model itself is incorrect.

restart;
with(plots):
Q := 1000: a := .9: b := 0.8e-3: mu := 0.247e-2: k := .2: y := 0.2e-1: e := .5: g := .5: T := 8: n := 100:
sigma[1] := 0.5e-1: sigma[2] := 0.9e-1: alpha[1] := 0.52e-2: alpha[2] := 0.52e-2: delta[1] := .8: delta[2] := .904: delta[3] := .8: c[1] := 50: c[2] := 250: c[3] := 50: w[1] := 140: w[2] := 130: w[3] := 150: w[4] := 160:
u[1] := min(max(0, z), 1): z := (a*i[p](t)*(i[p](t)+i[pm](t))*(lambda[4](t)-lambda[3](t))+a*s(t)*(i[p](t)+i[pm](t))*(lambda[1](t)-lambda[2](t)))/(n.w[1]): u[2] := min(max(0, C), 1): C := (b*(i[m](t)+i[pm](t))*s(t)*(lambda[3](t)-lambda[1](t))+b*(i[m](t)+i[pm](t))*i[p](t)*(lambda[4](t)-lambda[2](t)))/(n.w[2]): u[3] := min(max(0, j), 1): j := (i[p](t)*lambda[2](t)+i[pm](t)*(lambda[4](t)-lambda[7](t))-(i[p](t)+i[pm](t))*lambda[5](t))/w[3]: u[4] := min(max(0, o), 1): o := (i[m](t)*lambda[3](t)-i[pm](t)*(lambda[4](t)-lambda[7](t))-(i[m](t)+i[pm](t))*lambda[6](t))/w[4]: u[2] := 0: u[3] := 0:

sys := diff(s(t), t) = Q+delta[1]*r[p](t)+delta[2]*r[m](t)+delta[3]*r[pm](t)-(a*(1-u[1])*(i[p](t)+i[pm](t))/n+b*(1-u[2])*(i[m](t)+i[pm](t))/n+mu)*s(t), diff(i[p](t), t) = (1-u[1])*a*(i[p](t)+i[pm](t))*s(t)/n-(1-u[2])*b*(i[m](t)+i[pm](t))*i[p](t)/n-(sigma[1]+u[3])*i[p](t)-(alpha[1]+mu)*i[p](t), diff(i[m](t), t) = (1-u[2])*b*(i[m](t)+i[pm](t))*s(t)/n-(1-u[1])*a*(i[p](t)+i[pm](t))*i[m](t)/n-(sigma[2]+u[4])*i[m](t)-(alpha[2]+mu)*i[m](t), diff(i[pm](t), t) = (1-u[2])*b*(i[m](t)+i[pm](t))*i[p](t)/n+(1-u[1])*a*(i[p](t)+i[pm](t))*i[m](t)/n-(y+u[3]+u[4])*i[pm](t)-(alpha[1]+alpha[2]+mu)*i[pm](t), diff(r[p](t), t) = (sigma[1]+u[3])*i[p](t)+(e*y+u[3])*i[pm](t)-(delta[1]+mu)*r[p](t), diff(r[m](t), t) = (sigma[2]+u[4])*i[m](t)+(y*g*(1-e)+u[4])*i[pm](t)-(delta[2]+mu)*r[m](t), diff(r[pm](t), t) = (y*(1-g)*(1-e)+u[3]+u[4])*i[pm](t)-(delta[3]+mu)*r[pm](t), diff(lambda[1](t), t) = lambda[1](t)*(a*(1-u[1])*(i[p](t)+i[pm](t))/n+b*(1-u[2])*(i[m](t)+i[pm](t))/n+mu)-lambda[2](t)*(1-u[1])*a*(i[p](t)+i[pm](t))/n-lambda[3](t)*(1-u[2])*b*(i[m](t)+i[pm](t))/n, diff(lambda[2](t), t) = -c[1]+lambda[1](t)*(1-u[1])*a*s(t)/n-lambda[2](t)*((1-u[1])*a*s(t)/n+b*(1-u[2])*(i[m](t)+i[pm](t))/n+sigma[1]+u[3]+alpha[1]+mu)-lambda[4](t)*(b*(1-u[2])*(i[m](t)+i[pm](t))/n+(1-u[1])*a*i[m](t)/n)-lambda[5](t)*(sigma[1]+u[3]), diff(lambda[3](t), t) = -c[2]+lambda[1](t)*(1-u[2])*b*s(t)/n+lambda[2](t)*(1-u[2])*b*i[p](t)/n-lambda[3](t)*(u[4]+alpha[2]+sigma[2]+a*(1-u[1])*(i[p](t)+i[pm](t))/n-(1-u[2])*b*s(t)/n)-lambda[4](t)*((1-u[2])*b*i[p](t)/n+a*(1-u[1])*(i[p](t)+i[pm](t))/n)-lambda[6](t)*(sigma[2]+u[4]), diff(lambda[4](t), t) = -c[3]+lambda[1](t)*((1-u[1])*a*s(t)/n+(1-u[2])*b*s(t)/n)-lambda[2](t)*((1-u[2])*b*i[p](t)/n+(1-u[1])*a*s(t)/n)-lambda[2](t)*((1-u[2])*a*i[m](t)/n+b*(1-u[1])*s(t)/n)-lambda[4](t)*(y+u[3]+u[4]+alpha[1]+alpha[2]+mu-(1-u[2])*b*i[p](t)/n-b*(1-u[1])*i[m](t)/n)-lambda[5](t)*(e*y+u[3])-lambda[6](t)*(y*g*(1-e)+u[4])-lambda[7](t)*(y*(1-g)*(1-e)+u[3]+u[4]), diff(lambda[5](t), t) = -lambda[1](t)*delta[1]+lambda[5](t)*(delta[1]+mu), diff(lambda[6](t), t) = -lambda[1](t)*delta[2]+lambda[6](t)*(delta[2]+mu), diff(lambda[7](t), t) = -lambda[1](t)*delta[3]+lambda[7](t)*(delta[3]+mu), i[p](0) = 300, i[m](0) = 200, i[pm](0) = 150, r[p](0) = 200, r[m](0) = 150, r[pm](0) = 150, s(0) = 1000, lambda[1](T) = 0, lambda[2](T) = 0, lambda[3](T) = 0, lambda[4](T) = 0, lambda[5](T) = 0, lambda[6](T) = 0, lambda[7](T) = 0:
p1 := dsolve({sys}, type = numeric, method = bvp[midrich], abserr = 0.1e-3,maxmesh=2400);

Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging
 

Rewrite your equation as  k*x = 3/2*arctan(x)  and look at 2 graphs  y = k*x  and  y = 3/2*arctan(x) . Since the derivative of the second function at x=0 is 3/2, we get the obvious answer:   k∈(-∞,0]∪[3/2,+∞) . Below you can see the plot for k=3/2 :

 plot([3/2*x, 3/2*arctan(x)], x=-5..5, -5..5, color=[red,blue]);  

Edit. The solution above is true for the case b = 0. For arbitrary b<>0, we rewrite the equation  f(x)=b  as  k*x - b = 3/2*arctan(x) . Using the properties of functions and geomeric considerations, we get the final answer (k0(b) is the value of  k  for which for a given  b  the graphs touch each other ):

For  -3*Pi/4<b<0  or  0<b<3*Pi/4 :   k∈(-∞, 0] ∪ [k0(b), +∞) 
For  b<=-3*Pi/4  or  b>=3*Pi/4 :   k∈(-∞, 0) ∪ (0, +∞)

Below we can see the animation for the case b = -1. Note that the value of k0(b) can only be found numerically.

restart:
fsolve({k0*x+1=3*arctan(x)*(1/2),k0=3/2*1/(1+x^2)});
plots:-animate(plot,[[-3*Pi/4,3*Pi/4,3*arctan(x)*(1/2),[t*cos(s),1+t*sin(s),t=-5..5]], x=-5..5,-5..5,color=[black$2,blue,red], thickness=[1$2,2$2], linestyle=[3$2,1$2]], s=arctan(eval(k0,%))..Pi, frames=100);
                  {k0 = 0.3313385522, x = 1.878055290}
                  

See help on the  DEtools:-DEplot  command.

I think that it makes sense to write procedures for calculating the equivalent resistance for the general case, that is, for any number of resistances and their arbitrary connection. Below by combining  ser  and  par  procedures, we can calculate the equivalent resistance for more complex circuits:

Download circuits.mw

See help on the  GraphTheory:-CycleBasis  command.

The dimension of the solution of a system is equal to the number of arbitrary constants in its solution. But in this example, the solution depends on an arbitrary smooth function  x[1](t) . But the dimension of the space of such functions is equal to infinity:

restart;
Sys:={diff(x[1](t),t)+diff(x[2](t),t)+2*diff(x[3](t),t)=0, x[1](t)-2*x[2](t)-x[3](t)=0, 2*x[1](t)+x[2](t)+3*x[3](t)=0, -x[1](t)+2*x[2](t)+x[3](t)=0};
dsolve(Sys);

The  P  procedure solves the problem:

P:=(A::Matrix)->nops(op(2,A));

Example of use:

A:=<0,1,0,1; -1,3,0,0; 0,0,1,0>;
P(A);
                 

Here is another simple method using  plots:-animate  command. 40 frames for the ellipse (0<e<1, red color), 10 frames for the parabola (e=1, green color) and 40 frames for the hyperbola (1<e<=1.8, blue color) :

restart;
L:=[seq(n*0.025,n=0..39),1$10,seq(1.02+n*0.02,n=0..39)]:
plots:-animate(plot,[1/(1-e*cos(phi)), phi=0..2*Pi, color=piecewise(e<1,red,e=1,green,blue), thickness=2, coords=polar, scaling=constrained, view=[-10..15,-5..5], discont=true], e=L, size=[1000,300]);

restart;

expr1:=-lambda-(1/2)*kappa__c-gamma__p-(1/2)*sqrt(-16*N*g^2+4*lambda^2-8*lambda*gamma__p+4*lambda*kappa__c+4*gamma__p^2-4*gamma__p*kappa__c+kappa__c^2);

subsindets(expr1, sqrt, t->Student:-Precalculus:-CompleteSquare(t));
expand(%, op([4,2,1,1], %));

restart;
'eval'(diff(f(r),r), r=C)=0;

Edit.

The  irem  command does not work for symbolic variables  so we use the  isolve  command instead. The following code returns all pairs  [k,b]  for which a solution exists. Total list  L  of all the solutions contains 4000 solutions.

restart;
N:=0:
for k from 1 to 9999 do
sol:=isolve(k*b=m*10000+2391);
if sol<>NULL then N:=N+1; L[N]:=[k,eval(rhs(sol[1]),_Z1=0)]; fi;
od:
L:=convert(L, list);
nops(L);

Your pair  [297, 9503]  is also present in the list  L .

If you just need to animate a rolling disk, then the easiest way is to use the tools of  plottools  and  plots  packages.

Example:

restart;

F:=proc(t)
local Disk, Circle, Radius, Wheel, R;
uses plottools, plots;
Disk:=disk([0,1], color=yellow);
Circle:=circle([0,1], color=blue, thickness=3);
Radius:=line([0,0], [0,1], color=red, thickness=3);
Wheel:=display(Disk,Circle,Radius);
R:=rotate(Wheel, -t, [0,1]);
display(translate(R, t, 0), scaling=constrained);
end proc:

plots:-animate(F,[t], t=0..2*Pi, frames=120, size=[1000,400]);

We get the list of 5 solutions for this system. Using the  eval  command we can easily get the values of the expression  f  for each solution:

Download example_new.mw

Should be

Here is the plot of integration domain:

plot(s, s = 2 .. 4, filled, view = [0 .. 5, 0 .. 5], labels=[s,u]);

Download order_integration_new.mw

Addition. This integral is easy to calculate symbolically in two ways:

Download order_integration_new1.mw