The library must exist, or create a new one.

restart;
MyPackage := module() export f1, f2; local loc1; option package; 
  f1 := proc() loc1 end proc; 
  f2 := proc( v ) loc1 := v end proc; 
  loc1 := 2; 
end module:
mylib:="C:/tmp/test.mla";
LibraryTools:-Create(mylib);

#                   mylib := "C:/tmp/test.mla"

LibraryTools:-Save(MyPackage, mylib);
LibraryTools:-ShowContents(mylib); # optional
#    [["MyPackage.m", [2021, 9, 15, 12, 20, 8], 41984, 126]]

restart;
libname:=libname, "C:/tmp/test.mla":
with(MyPackage);
#                            [f1, f2]

f1();
#                               2

f2(123);
#                              123

f1();
#                              123

restart;
nat_param:=proc(ff, xx ::(name=anything), yy::(name=anything), t::name)
  local K,xt,yt, ode;
  xt:= diff(ff, op(1,xx)); 
  yt:= diff(ff, op(1,yy));
  K := 1/sqrt(xt^2 + yt^2);
  ode := {'diff'( op(1,xx),t ) = -K*yt, 'diff'( op(1,yy),t ) = K*xt};
  ode := subs( [op(1,xx)=op(1,xx)(t), op(1,yy)=op(1,yy)(t) ], ode);
  dsolve(ode union {op(1,xx)(0)=op(2,xx), op(1,yy)(0)=op(2,yy)},numeric, output=listprocedure);
  eval([lhs(xx(t)),lhs(yy(t))],%)[]
end:

f:=x1^4 + x2^4 + 0.4*sin(7*x1) + 0.3*sin(4*Pi*x2) - 1: 
xx,yy:=nat_param(f, x1=0,x2=1, t):
L:=fsolve(xx(t)=xx(0),t=8..10): N:=6:
plots:-display(seq(plot([xx,yy,k*L/N..(k+1)*L/N],color=COLOR(HUE,k/(N+2))), k=0..N-1), thickness=6);

op(1, xLp_tmp)

It's a bug (but actually the generators are correct).

restart;
with(GroupTheory):
H:=PermutationGroup({[[2, 3], [4, 5]], [[2, 5], [3, 4]]}, degree = 5, supergroup = PermutationGroup({[[2, 3, 4]], [[1, 2], [4, 5]]}, degree = 5));
#                    H := ((23)(45)(25)(34))

MinPermRepDegree(H);
#"GENS" = [[Perm([[1, 2]])], [Perm([[3, 4]])]]

#Error, (in GroupTheory:-Subgroup) not all the provided generators belong to the group

Workaround:

GroupOrder(H);
 #                              4

H1:=DirectProduct(CyclicGroup(2)$2):
AreIsomorphic(H, H1);
#                              true

MinimumPermutationRepresentationDegree(H1);
#"GENS" = [[Perm([[1, 2]])], [Perm([[3, 4]])]]

#                               4

The command GroupTheory:-IdentifySmallGroup identifies a group returning (n,d) where n is its order and d an integer (1 <= d <= NumGroups(n)). Then inspect https://en.wikipedia.org/wiki/List_of_small_groups

The computation is slow on my (old) computer. It needs several hours.
To  interrupt and then continue later, just save the partial results. 

restart;
read "d:/temp/Jfile.mpl";
max(indices(J,nolist));
smax := 20;  R := 1680;  k := 1/4;  gam := 1/2;    vmax:=R/2;
n:=0;
Digits:=15:
f:=proc(v,s) 
  evalf(Re(  sinh(s)*coth(1/2*s)^(2*I/k)/(gam-I*k*cosh(s))^5*exp(-(k^2*sinh(s)^2/(gam-I*
  k*cosh(s))+I*k*(1-cosh(s)))*v)*Hypergeom([I/k],[1],I*k*v)*v^n  ))
end:
V0:=0;  dV:=20;  # <-----
for V from V0 to vmax-dV by dV do
  J[V]:=evalf(1/R * Int(f, [V..V+dV, 0..smax], method=_CubaCuhre, epsilon=1e-5));
od;
save J, "d:/temp/Jfile.mpl";
add(entries(J,nolist));
max(indices(J,nolist));

Similar for the imaginary part if needed.

Download SimpJP-vv.mw

F = 0, G = 0 is obviously a solution (gamma = anything).
Depending on your coefficients, this could be the only one.

A:=<-2,3,-5>: B:=<-6,1,-1>: C:=<2,-3,7>:
with(LinearAlgebra): local D:
k:=Norm(B-A,2) / Norm(C-A,2):
D:=(B + k*C)/(1+k);
eqAD:=A + t*(D-A); # parametric equation

The first problem is that you cannot have t=1:

fN(0.5, 1);
Error, (in solnproc) unable to compute solution for t>HFloat(0.2):
Newton iteration is not converging

So, let's take  t := 0.1;

It is not possible to compute the indefinite integral numerically. In your case, what you want is a definite integral.

t := 0.1;
                            t := 0.1

A1:=x*R(z)*R(z)*(fN)(x, t);
                A1 := 0.8692871388 x fN(x, 0.1)

A2:= X -> int(A1, x=0..X):
A2(0); # 0, as you want
                               0.

plot(A2, 0..1);

Download Worksheet_1_vv.mw

If you really want to create the (global) m||i variables even when they exist, just replace

M := [seq(m || i, i = 2 .. l)];

with

M := [seq(evaln(m || i), i = 2 .. l)];

It's a quadratic equation. Just use solve.

solve((a*y-x*(a+1))/(x^(2)-b*a^(2)*(y-x)^(2))=A, y);