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.

@Kitonum I don't think that this is a bug. A named procedure evaluates to its name.
So what you need to do is
 

is(eval(f1)=eval(f2));

which returns true.

Why do you want Maple 6? Since then there have been 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 2015, 2016, and 2017.

In any case the correct address for your inquiry is the maplesoft website:
maplesoft.com

I ran your code in Maple 2017.3 and didn't have any problem.
The plot:

What is the value of N?
The variable tot is not really used for anything, just assigned to 0.
Did you mean this code to run immediately after a restart? 

@Mariusz Iwaniuk Yes, only minute changes need to be made to the answer I gave in your link:
 

restart;
ode:=diff(y(x), x$3)+diff(y(x), x$2)+y(x)=1; 
#Now only 2 points, but with an additional requirement at 0, in this case D(y)(0) = a, where a is to be determined so that y(1) = 0:
#
p:=a->dsolve({ode,y(0)=0,y(2)=1,D(y)(0)=a},numeric,output=listprocedure);
#
p1:=proc(a) local Y1; Y1:=subs(p(a),y(x)); Y1(1) end proc;
#Just exploring the ground:
plots:-odeplot(p(-1),[x,y(x)],0..2);
#Now looking for a zero for p1, i.e. a value for a for which y(1) = 0 using a = -1 as an initial guess:
A:=fsolve(p1,-1);
p(A)(1);
plots:-odeplot(p(A),[x,y(x)],0..2);