Just to differentiate the results for different n. For the n 's I tested the curves of the original problem are the same, so I added "something" to make them different - as a test. Completely arbitrary! You might as well add any other function of n to see that differerent n's give different results. When you then eliminate that additional term I think you will see only one line, indicating that you n has a neglibible influence on your result.

Perhaps you meant something like this

pp[nn_] := Module[{},
  n = nn;
  y = Sum[x^i*a[i], {i, 2, n}];
  th = D[y, x];
  M = EI*D[y, {x, 2}];
  V = EI*D[y, {x, 3}];
  PE = EI/2*Integrate[D[y, {x, 2}]^2, {x, 0, L}] + P*(y /. x -> L);
  Eq = Table[D[PE, a[i]], {i, 1, n}];
  sol = Solve[Eq == 0, Array[a, n]];
  FullSimplify[M /. sol[[1]]];
  FullSimplify[V /. sol[[1]]];
  EI = 100;
  P = 1;
  L = 10;
  Plot[{M /. sol, -100, 100}, {x, 0, L}, 
   AxesLabel -> {"x (m)", "M (N.m.)"}]
  ]

Off[Solve::"svars"]
Show[Table[pp[j], {j, 3, 6}], PlotRange -> {-15, 15}]

It seems that all lines coincide.

If you change (just as a test)

Plot[{M /. sol, -100, 100}, {x, 0, L}, AxesLabel -> {"x (m)", "M (N.m.)"}]

to

Plot[{ n/2 +M /. sol, -100, 100}, {x, 0, L}, AxesLabel -> {"x (m)", "M (N.m.)"}]

you will see differences

Crossposted here.