Torsten is right. There is no way you will get a solution with a large step size (h=5). I checked out your code and it works fine for mu=1 with a step size of h=0.1
I would recommend the following two changes;
1) add the command
at the start of the script.
2) I would also recommending pre-allocating the yb2 vector. you can do this by changing
this will make yb2 and yb the same size, and they will both have a 1 in the first position.
Now, as you increase the value of mu, the system becomes more stiff. For example, using a step size of h=0.1 I can get a solution for mu =2 and mu=3, but at mu=4 the function blows up after about 14 seconds. (it is numerically unstable). Now you have to reduce the step size to improve the calculation stability. In order to get a solution with mu=1000, and to execute for 3000 seconds without blowing up, the step size must be extremely small (10^-6 or even smaller). My problem is that this is taxing the capabilities of my computer and I cannot find the solution.
A major part of the problem is the use of the Euler integration method. This method is not well suited for stiff systems (perhaps that is what this lesson is about). For this problem, I would switch to a Runge-Kutta integrator and try my favorite [5,6] Runge Kutta method. This would allow the use of a much larger step size (probably 3 orders of magnitude) compared to what is required for Euler integration.