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

 Maple simply shows the first few and last few digits of the numerator and denominator of this gigantic fraction. If you are confused by the incorrect seventh digit in the numerical calculation, then simply increase Digits . It is never guaranteed that when calculating (with  Digits=10 by default) all 10 digits are correct.
 

If you raise the number close to 1 to a very high degree with the accuracy by default, an even worse situation may occur (the correct answer is close to the number e);

restart;
(1+1./10^10)^(10^10);
                                           1.000000000


Edit.

Download arithmetic_new.mw

restart;

f:=(x,y)->x*y^2+1:  # Specify the function
TP:=f(2,1)+D[1](f)(2,1)*(x-2)+D[2](f)(2,1)*(y-1);  # Find the tangent plane at the point (2,1)

plot3d([f(x,y), TP], x=2..3, y=1..3, color=[yellow,blue], axes=normal);

plot3d([f(x,y), TP], x=2..2.1, y=1..1.2, color=[yellow,blue], axes=normal);

plot3d([f(x,y), TP], x=2..2.01, y=1..1.02, color=[yellow,blue], axes=normal);

Use the  solve  command for this. We get 2 solutions:

solve({x+y=1, y=x^2-5});