@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);

@vv Yes, but that is to be expected since [0.3333] will remain [0.3333] after evaluation, whereas the evaluation of [0.3333]+[0] will involve doing 0.3333 + 0, which (with Digits=3) will result in 0.333.

@vv I see! Thank you.
It is important for this example that u is a name assigned to 0.3333 (when Digits=3).
When you take the number itself:

[0.3333+0] = [0.3333] + [0];

you get 0.333 = 0.333.

To make automatic simplification come first on both sides in the name assigned case we can do
 

'[u+0]' = '[u] + [0]';

which evaluates to [u] = [u]  and after %; to .3333 = .3333.
 

@vv Sorry, I don't get it.
Can you shows us such an example?

I'm not able to explain everything here, but will just mention two points that play a role:
1. The input of a floating point number alone as in 1.23456;  leaves all decimals even though Digits < 6. Only if arithmetic with floats come into play as in 1.23456 + 0; is Digits respected.
2. Automatic simplification kicks in here before arithmetic with floats:
 

1.23456;
%+0;

and here:
 

u:=1.23456;
u+0;

Thus the output as a consequence of point 1 is 1.23456 in both cases.