Hi Mohga,
Your result does not match the expected result because you didn't graph the same values for Z. You've also got a vulnerability by using machine numbers in your table, which will often cause Mathematica to lose internal precision because it will downgrade an expression result to account for the numerical properties of machine precision values. Try the following, which is a minimal edit of your version taking the above into account:
myprecision = 30;
Clear[y, x];
Clear[Derivative];
soln = ParametricNDSolve[{y''[
x] == (1 - z) Sin[
y[x]] + ((z*y[x])/(2 \[Pi]^2)) (y[x] -
2 \[Pi]) (y[x] - \[Pi]), y[15] == 2 Pi, y[-15] == 0},
y, {x, -15, 15}, {z}, Method -> "BDF", MaxSteps -> Infinity, WorkingPrecision -> myprecision];
Plot[Evaluate[
Table[y[SetPrecision[z, myprecision]][SetPrecision[x, myprecision]]/
Pi /. soln, {z, {0, 0.5, 1, 2, 2.5, 2.609}}]], {x, -15, 15},
PlotRange -> All]
You can generally increase performance substantially by tuning AccuracyGoal, PrecisionGoal, and the step/recursion parameters for your chosen method. I usually start with conservative values, check that the result behaves well as I vary the values, then tune for performance for my particular situation.