@Adam Ledger If you do this:
 

subs(_Z=myZ,RootOf(_Z^2-2));
%;

you will see that you are back where you started.
Why mess with the design of RootOf?

@mmcdara A way out seems to be to forget about the file TESTproduce.mw and use this version of TEST_Examples.mw instead of the original. No use of unwith, but the reading and modification of the .mla file takes place here:
 

restart;
with(TEST);
TestProcedure();
####### Now change the mpl text file to version 2 and then do:
####### No restart (the point).
read "C:/Users/Bruger/maple/toolbox/TEST/TEST.mpl";
savelibname:="C:/Users/Bruger/maple/toolbox/TEST/lib";
savelib(TEST);
#######
TestProcedure();

@mmcdara I have uploaded 3 files.
TEST.mpl is a text file, which happens to be created by Export as Maple Input (.mpl).
TESTproduce.mw is a worksheet which reads TEST.mpl and produces content for TEST.mla.
TEST_Examples.mw is a worksheet with one restart only.
The assumtion is that the archive TEST.mla has been created in the place indicated.
###
At start I used TEST.mpl as written (Version 1).
In Testproduce.mw I check that it works.
In TEST_Examples.mw I start from the beginning and continue to the line " ####### Now change ... "
Then I change TEST.mpl to Version 2.
Go to Testproduce.mw and check that it works: now printing "Version 2".
Finally go to TEST_Examples.mw (without restarting: the point of the exercise).
Start with unwith(TEST); and continue.
You get Version 1, which didn't surprise me at all.

TESTproduce.mw
TEST_Examples.mw
MaplePrimes won't let you upload .mpl files, but TEST.mpl is just this and no more:

TEST:=module() option package; export TestProcedure;
  TestProcedure:=proc() print("Version 1"); 123456789 end proc
end module;

@Rouben Rostamian  _Z is used by RootOf in e.g. this explicit call to RootOf:
 

RootOf(x^2-2,x);

RootOf( _Z^2 - 2) is returned and _Z is the global _Z.
Compare this to the following

restart;
sol1:=solve(sin(x)=0,allsolutions);
sol2:=solve(sin(x)=0,allsolutions);
S:=indets({sol1,sol2},And(name,Not(constant)));
type~(S,`local`); #{true}
_Z1:=34; #No problem _Z1 is the global _Z1
S; # No change

So it is good that the help page for RootOf points out that the _Z appearing there is the global _Z.
After all _Z could have been local:
 

restart;
R:=proc(pol,x) local _Z; rootof(subs(x=_Z,pol)) end proc;
R(x^2-2,x);

So the question is rather: why is it good (or even necessary) for _Z to be global in RootOf?

@Rouben Rostamian  You wrote:
" The error message that you have shown does not complain about "global".  It complains about "protected".  It may be a global variable, but that's beside the point. "

To that I can only repeat that the name _Z in RootOf is indeed a global variable (and it is also protected). Why protect a local variable? Since the question raised by Adam Ledger was about _Z being global, I cannot see that it is besides the point (as you say).

You may also try
 

res:=solve(x^5+5*x^4+4*x^3+2*x^2+x+9);
indets(op(1,res[1]));
type(op(%),`global`); #Answer true

As is well-known dsolve uses _C1, _C2, ... for arbitrary constants. They are globals, but not protected.
 

sol:=dsolve(diff(y(x),x)=y(x));
S:=indets(sol,name) minus {x};
type(op(S),`global`); # true
type(op(S),protected); #false
_C1:=89;
dsolve(diff(y(x),x)=y(x)); #Now using _C2 because of the assignment made above!

Certainly unwith works with user generated packages. I just tried. No problem.
When I make changes to a procedure in a package (happens often and is done in a text editor) then I regenerate the mla file in a worksheet made exclusively for that purpose.
In my Examples worksheet I begin with a restart.
All works well.
So your situation is that you don't want to begin with a restart in your "Examples" worksheet.
Is this close enough to your situation?
####
You want to do this in sequence as written
 

unwith(MyModule):
with(MyModule);   

and expect the first MyModule to refer to the old version and the second MyModule to refer to the new version. Correct? 

@Rouben Rostamian  Well if _Z wasn't global I wouldn't get the complaint I get here:
 

res:=solve(x^5+5*x^4+4*x^3+2*x^2+x+9);
_Z:=987;

We get the error:

Error, attempting to assign to `_Z` which is protected.  Try declaring `local _Z`; see ?protect for details.
We can continue with (but I wouldn't recommend that)
 

unprotect(_Z);
_Z:=987;

Having said all that I agree with everything else you said in your answer.

@student_md There is no way of solving the system numerically if F is not a known function or if F(t) doesn't evaluate to a number when given a number t, not in Maple nor in any other numerical program.

@tomleslie I don't know if Rouben Rostamian would like to be mistaken for me :-)

@student_md If alpha=0 then the system is very simple and can easily be solved analytically:

## The general solution in the alpha=0 case:
dsolve(eval(SYS,alpha=0));

Although you have a linear combination of sines and cosines (each periodic) the solution will not be periodic because the ratios of the periods are irrational.
For alpha small the behavior will look similar.
You could try the numerical approach given above with alpha = 0.01 as the only change (i.e., all other parameters=1).

@student_md Without using Maple your only way out is to use another program with the same capabilities.
I don't think you will be able to get an analytical solution by hand or in Maple or any other program.
Try the following in Maple, where I use the worksheet I uploaded:

## SYS known from worksheet
## paramvalues known from worksheet
## Isolating the highest derivatives (2nd order):
solve(SYS,diff~({V1,V2,V3,V4}(t),t,t));
## Picking the right hand side of the equation for V1''(t):
RHS1:=subs(%,diff(V1(t),t,t));
## Setting V2=V3=V4=0 gives us an ode in just V1.
## I'm also setting all parameters equal to 1:
ode:=eval(diff(V1(t),t,t)=eval(RHS1,{V2=0,V3=0,V4=0}),paramvalues);
## Presumably ode should be easier to solve than the whole system.
## But try even without initial conditions
dsolve(ode);

The answer is basically a restatement of the problem in the form of a DESol structure:

There is a lot of numerical code available in different languages, so you have a lot of choice in solving numerically.
In Maple the default method is Runge-Kutta-Fehlberg (rkf45). Other methods are available.
See the help page ?dsolve,numeric,ivp where also other options are mentioned.

@rahinui You should really be using pdsolve. That is made for pdes.
In this case it doesn't help, though.

@Kitonum Briefly, I had a reply reporting a length of 3832 for EXPR2.
I was using a beta version.
In Maple 2017.3 I get 3766.
I shall report that beta behavior.
 

@Haile With u[2] = 0 use abserr = 1e-4  (also written 0.1e-3):
 

p1:=dsolve({sys}, numeric, abserr = 0.1e-3,continuation=cont,maxmesh=2048);

@tsunamiBTP Your last definition of T1 suffers from the problem I described earlier.
You have   T1 := m-> 2*sin(alpha)*(diff(S11, t));
When T1 is applied as in
T1(7);
the first that happens is that the input (7) is evaluated (there is not much to do, 7 is 7 after all).
After that the formal parameter m is replaced by 7 in the unevaluated body of the procedure, i.e. in 2*sin(alpha)*(diff(S11, t)).
But in fact there is no m there to be seen, so no change is made.
You need unapply!
##
When should we consider two procedures equal?
My answer would be that they should have exactly the same content. Yours don't.
Consider the following three simple procedures:

p1:=proc(x) local y; y:=x^2; y end proc;
p2:=proc(x) x^2 end proc;
is(eval(p1)=eval(p2)); # false
p3:=proc(y) y^2 end proc; # p1 with just a name change
is(eval(p3)=eval(p2)); # true

Let me add that I'm no big fan of using is and would never in my wildest dream have thought of applying is to compare procedures.