It works, but the results must be wrong, because when I choose Sigma to be small and increase Mu drastically, I should get something which looks like an exponentially increasing function, because the function is driven heavily by the drift. However, according to my code this is not the case. So something must be wrong.

Does anybody have a solution when I want to manipulate the start value (S1), Sigma, and Mu, for one and for multiple paths?

Just generic remarks : I don't think your solution for S(t) is correct. You don't accumulate the "random stuff" only and then add the drift bit at the end. Also, the section regarding multiple paths doesn't seem relevant at this point.
–
b.gatessucksMay 24 '13 at 19:49

@Milan Ivica I really don't understand why you keep insisting in developing your own proprietary code (or use someone's code), when Mathematica already have all this GBM stuff internally... I think in my answer I've already shown exactly what you want for the 1-case and multiple-case paths...
–
RodMay 24 '13 at 20:24

Mathematica is a registered trademark of Wolfram Research, Inc. While the mark is used herein with the limited permission of Wolfram Research, Stack Exchange and this site disclaim all affiliation therewith.