In principle, Student:-Precalculus:-CompleteSquare should work but in this case it complicates things.
I prefer to use my generalized version:

SQR:=proc(P::polynom(anything,x), x::name)
local n:=degree(P,x)/2, q,r,Q,R,k;
if not(type(n,posint)) then error "degree(P) must be even" fi;
Q:=add(q[k]*x^k,k=0..n-1) + sqrt(lcoeff(P,x))*x^n;
R:=add(r[k]*x^k,k=0..n-1); 
solve({coeffs(expand(P-Q^2-R),x)}, {seq(q[k],k=0..n-1),seq(r[k],k=0..n-1)});
eval(Q,%)^2+eval(R,%)
end:

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

evalindets(expr1, sqrt, t -> SQR( op(1,t), lambda)^op(2,t));

{seq(msolve({a = k, a*b = 2391}, 10000), k = 1 .. 9999)};

with(ImageTools): 
img0:=Read("smile.jpg"):   str:="Just a test!":
img1:=Scale(img0,1..50):
img:=Create(400,400,channels=3,background=white):
simg:=(x,y) -> SetSubImage(img,img1, min(max(floor(400-400*y-25),1),350), min(max(floor(400*x-25),1),350)):
Explore( plot(x^2,x=0..1, background=simg(a,a^2), axes=none, title=str[1..floor(15*a)]), a=0..1., animate,loop,autorun);

(The file "smile.jpg"  must be in the current directory).

Your system has a sigularity at t=0.
Just replace in ICS:  S(0) = 30 with say S(0.001) = 30  etc. It will work.
 

It was answered recently here.

You cannot use the operator * as a name. So, use  beta[`*`], or even better `beta__*`.

Also, remove the last comma in params.
 

In Maple there are many options for numerical integration (including compiled external libraries) which usually can increase considerably the speed. You should choose the proper ones for your case. Read the help pages:

?evalf,int

and also the provided links.

# The coordonates for the center of B is (x, r+d).
x^2 + (d+r)^2 = (R+r)^2;
solve(%,x);
u := factor(%[1]);

# The coordinates of the contact are:
R/(R+r) * <u, d+r>;

a:=1; b:=2; c:=3;  # semi axes
plot3d([a*cos(theta)*sin(phi),b*sin(theta)*sin(phi),c*cos(phi)],
theta=0..2*Pi,phi=0..Pi,scaling=constrained,axes=boxed);

For a=b=c=1 (as in your code) you have a sphere.

restart;   # Continuous roots w1(k)..w5(k) of a polynomial (wrt parameter k in 0..1)
with(plots):
f := -135/4*w^5+369/16*w^3*k^2+47/4*I*w^4-93/16*I*w^2*k^2+w^3-2/3*w^3*k*B-27/16*k^4*w+3/16*I*k^4-1/3*w*k^2+2/9*k^3*w*B:
B0:=1:
ini:=fsolve(eval(f,[B=B0,k=1/2]),w):
ode:=diff(eval(f,[w=W(k),B=B0]),k):
solode:=seq(dsolve({ode,W(1/2)=ii},numeric),ii=[ini]):
display(Matrix(2,5,
        [ [seq( odeplot(solode[i], [k,Re(W(k))], k=0..1,title=w[i]), i=1..5)],
          [seq( odeplot(solode[i], [k,Im(W(k))], k=0..1), i=1..5)] ]
),thickness=3);

Download PolynomialFit-simplex.mw

I seems that a polynomial of degree <=13 satisfying the "marked" conditions in the provided graph and having also a similar "shape" (i.e. same sign for the derivative) does not exist.

N:=13:
t:= unapply(add(x^k*u[k],k=0..N),x):
sol:=solve([
        t(-7/2)=-5,
        t(-2)  = 2, 
        t(0)   = -2, 
        t(1)   = 0, 
        t(5/2) = -3, 
        t(4)=3, 
        D(t)(-2)=0,
        D(t)(0)=0,
        D(t)(1)=0,
        D(t)(5/2)=0,
        D(t)(-7/2)=a,  # a>0
        D(t)(4)=b,     # b>0
        t(-1)=c,       # -1<c<1
        t(7/2)=d       # -1<d<1
      ], [seq(u[i],i=0..N)] 
     
)[]:
P:=eval(t(x),sol):
#indets(P);
Explore(plot(P, x=-3.5..4), parameters=[[a=0. .. 100],[b=0. .. 100], [c=-1. .. 1],[d=-1. .. 1]]);

residue connot compute the residue at a pole with symbolic order.
It's easy to write a procedure in this case provided that you know the order of the pole. It is based on Maple's ability to compute the derivative of a symbolic order, but it will fail for complex expressions.

Res:=(f,z,a,p) -> limit(diff(f*(z-a)^p,[z$(p-1)]), z=a)/(p-1)!:
# p = the order of the pole z=a

f := exp(-z/A)/z^(2 + n):
Res(f,z,0,n+2);
        

Try to increase Digits. For example,

restart;
Digits:=10;
n:=10;
x:=ln(1+0.1^n);
x^x;

==> Float(undefined). Increasing Digits ==> OK.