You must use

Explore('a'[i], i = 1 .. 3);

otherwise accessing a[i]  with symbolic i, produces an error.

(Explore has special evaluation rules, the first argument is passed unevaluated).

pde:={diff(u(x,t),t) =K*diff(u(x,t),x, x), D[1](u)(0,t)=0}: 
# u(x,infinity) = infinity,  diff(u(x,t),x) =0 when x=0 for all t 
sol:=pdsolve(pde, generalsolution);

# ==> (changing the constants)

u(x,t) = A*( exp(c*(K*c*t+x))+exp(c*(K*c*t-x)) );

# A>0, c>0 [K>0]
# pdsolve says it's the general solution, but with pdsolve you can never know

SolveTools:-SemiAlgebraic([x^2+y^2+z^2-4*x+6*y-2*z-11 = 0, 2*x+2*y-z = -18], [x,y,z]);

       [[x = -4/3, y = -19/3, z = 8/3]]
              

You should post text, or preferably Maple code instead of low resolution pictures.

Anyway, you do not need Change if your f is smooth:
int(diff(f(s),s), s = a..b, continuous);
      - f(a) + f(b)

Chebyshev already had the smart idea in 1853, see here.

After a change of variables ( x+1 --> x), it's a binomial integral

Int( x^m * (A + B*x^n)^p, x);

where m=-3/2, n=3*b, p=-1/2

Using Chebyshev's result, assuming that b is rational (and a<>0, a<>1 of course), the integral is an elementary function iff
b = 1/(6*k), or b = 1/(6*k + 3),  where k is an integer.

Yes, the computation are very slow in the general case.

You can speed up considerably the computations if you use hardware floats.

If A is the 0-1 matrix, define B:=Matrix(A, datatype=float)   and compute C := B^k.
The result will be exact provided that the entries of C are < 2^53 =~ 9.0*10^15.

Another posibility is to use LinearAlgebra:-Modular.

Here I have assumed the usual multiplication. For a custom multiplication you may use LinearAlgebra:-Generic.

The level sets BS_Price = B (B constant) are 4-dimensional manifolds.
The equation can be solved e.g. wrt K. (It is convenient to use IERF the inverse erf function, see the code below).

To visualize the 3d sections, one can use Explore.

restart;
BS_Price=exp(-r*T)*(-(1/2)*erf((1/4)*sqrt(2)*(sigma^2*T-2*r*T+2*ln(K)-2*ln(S__0))/(sigma*sqrt(T)))+1/2);
IERF(2*B*exp(r*T)-1) = (-(1/2)*sigma^2*T+r*T-ln(K)+ln(S__0))*sqrt(2)/(2*sigma*sqrt(T));
EQK:=solve(%, K);
IERF:= z -> RootOf(erf(_Z)-z):
Explore(plot3d(EQK, r = 0 .. -ln(B)/T, T=0.5 .. 2, labels=[r,T,K]), 
parameters=[sigma=0.01 ..1, S__0=0.01 ..2, B=0.01  ..0.95]);

AFAIK combinat[randperm] exists in Maple versions much older than 13.

Maple can use theoretically up to 3.8*10^10 digits - if you have enough memory :-).
Your number has 83196715617456  digits (i.e. > 8*10^13) , so it will "overflow".

Using Maple + some elementary math yo can compute

a) the first n digits
b) the last n digits
for any reasonable n.
Try it.

A much more dificult task (probably impossible) would be to compute e.g. the 2 digits in the middle for this number.
If you can do it, please let us know!

1.  Use \\  or  /  as  directory separator.

2. The assignment operator is := 
so, use its := 1;

simplify((1/2)!);
  

You cannot simplify to 0 because it is not 0. But it is a constant.
You can check this in Maple, computing

diff(expr, x);   # 0

where expr is your expression. Do not forget to use u(x) and v(x)  instead of u, v.
 

Tunnell:=proc(N::posint)
local A, a,t, n, isq;
n:=mul(map(`^`~ @ `op`, map(t->[t[1],irem(t[2],2)], ifactors(N)[2])));
isq:=proc(n) local r:=isqrt(trunc(n)); `if`(r^2>n,r-1,r) end;

A:=proc(a,b,c)
  local x,y,z, r:=0;
  for z from 0 to isq(n/c) do
  for y from 0 to isq((n - c*z^2)/b) do
    x:=(n-c*z^2-b*y^2)/a;
    if issqr(x) then r:=r+`if`(z=0,1,2)*`if`(y=0,1,2)*`if`(x=0,1,2) fi
  od od: 
r 
end;

if    n::odd  and 2*A(1,2,32)=A(1,2,8) 
   or n::even and 2*A(2,8,64)=A(2,8,16)
then true else false fi
end:

select(Tunnell, [seq(1 .. 200)]);

[5, 6, 7, 13, 14, 15, 20, 21, 22, 23, 24, 28, 29, 30, 31, 34, 37, 38, 39, 41, 45, 46, 47, 52, 53, 54, 55, 56, 60, 61, 62, 63, 65, 69, 70, 71, 77, 78, 79, 80, 84, 85, 86, 87, 88, 92, 93, 94, 95, 96, 101, 102, 103, 109, 110, 111, 112, 116, 117, 118, 119, 120, 124, 125, 126, 127, 133, 134, 135, 136, 137, 138, 141, 142, 143, 145, 148, 149, 150, 151, 152, 154, 156, 157, 158, 159, 161, 164, 165, 166, 167, 173, 174, 175, 180, 181, 182, 183, 184, 188, 189, 190, 191, 194, 197, 198, 199]

Edit. The proc was optimized for speed.

A module is very convenient and has advantages, but it is not mandatory here.

with(Statistics):

Levy:= (m,s) -> Distribution(
                      PDF=(x->piecewise(x<0, 0, sqrt(s/2/Pi)*exp(-s/2/(x-m))/(x-m)^(3/2))), 
                      CDF= (x->erfc(sqrt(s/2/(x-m)))), 
                      RandomSample = (N->m +~ s /~ ( Quantile~(Normal(0, 1), Sample(Uniform(1/2, 1), N)) )^~2)
                    ):

X := RandomVariable(Levy(m,s)):
PDF(X, x);
Sample(X, 5);
Y := RandomVariable(Levy(1,2)):
PDF(Y, x);
Sample(Y, 6);