@Bendesarts Use output=listprocedure, and plot as a parametric curve:

sol:=dsolve({sys,Cinit},numeric,output=listprocedure):
Sol := subs( sol, [x(t),z(t)] ): X,Z := Sol[]:
plot( [X(t),Z(t),t=0..10],numpoints=200,color=blue,legend="z(x)");

@sunit Why Newton-Raphson and not just fsolve? What did you try? Submit a new question if your attempt doesn't work properly.

@sunit Why Newton-Raphson and not just fsolve? What did you try? Submit a new question if your attempt doesn't work properly.

@sunit Notice that I start with omega>0 to prevent division by zero:

pnts := [seq( [omega, fsolve(f)], omega=0.001..10,0.1 )]:
y := (omega,x) -> [x,x^2/omega]:
ypts := (y@op)~(pnts):
plot(ypts);

Explanation: I make y a function of two variables, returning a list consisting of the second variable and the calculated y-value.
Now you can not apply y directly to all elements of the list pnts by y~(pnts), because these elements are lists, so we make a composition of y and op, which extract the elements of the sublists of pnts.

@sunit Notice that I start with omega>0 to prevent division by zero:

pnts := [seq( [omega, fsolve(f)], omega=0.001..10,0.1 )]:
y := (omega,x) -> [x,x^2/omega]:
ypts := (y@op)~(pnts):
plot(ypts);

Explanation: I make y a function of two variables, returning a list consisting of the second variable and the calculated y-value.
Now you can not apply y directly to all elements of the list pnts by y~(pnts), because these elements are lists, so we make a composition of y and op, which extract the elements of the sublists of pnts.

@Alejandro Jakubi and @Carl Love 

This is not the first time that you teach me something I didn't know! I was aware of the distinction between x^(1/2) and x^.5. A demonstration:

Download Sqrt.mw

@Alejandro Jakubi and @Carl Love 

This is not the first time that you teach me something I didn't know! I was aware of the distinction between x^(1/2) and x^.5. A demonstration:

Download Sqrt.mw

@N00bstyle I added the command

print(subs(lengte = _lengte, eigenfrequentie));

in your procedure reductiefactorvoetgangerslengte. Now I get

lengte := 22;
reductiefactorvoetgangerslengte(lengte);
                0.0077 √ 234256     
Error, (in reductiefactorvoetgangerslengte) cannot determine if this expression is true or false: 0.769201511143595e-2*234256^(1/2) <= 1.25

but

evalb(0.77e-2*sqrt(234256)<=1.25);
                             false

@N00bstyle I added the command

print(subs(lengte = _lengte, eigenfrequentie));

in your procedure reductiefactorvoetgangerslengte. Now I get

lengte := 22;
reductiefactorvoetgangerslengte(lengte);
                0.0077 √ 234256     
Error, (in reductiefactorvoetgangerslengte) cannot determine if this expression is true or false: 0.769201511143595e-2*234256^(1/2) <= 1.25

but

evalb(0.77e-2*sqrt(234256)<=1.25);
                             false

... because f(1) = a+b, so perfectly well defined. But the right derivative only exists if a+b=1.

... because f(1) = a+b, so perfectly well defined. But the right derivative only exists if a+b=1.

Of course. I only wanted to show an inexpected pecularity!

My answer focuses on the question why youre original code didn't work. This results in a rather obscure program. This can be made more transparent whwn you realize that eq1 and eq2 are intended to be dependent on the variable a, so make these functions of a:

restart; B := 53*Pi*(1/180): e := 3:
eq1 := a -> -2.005689708*a:
eq2 := a -> -2.005689708*a+5.369606170:
f := proc (a)
  if a <= 0 then 0 
  elif a <= evalf(e*sin((1/2)*B)) then eq1(a)
  elif a <= evalf(2*e*sin((1/2)*B)) then eq2(a) 
  else 0 
  end if 
end proc:

My answer focuses on the question why youre original code didn't work. This results in a rather obscure program. This can be made more transparent whwn you realize that eq1 and eq2 are intended to be dependent on the variable a, so make these functions of a:

restart; B := 53*Pi*(1/180): e := 3:
eq1 := a -> -2.005689708*a:
eq2 := a -> -2.005689708*a+5.369606170:
f := proc (a)
  if a <= 0 then 0 
  elif a <= evalf(e*sin((1/2)*B)) then eq1(a)
  elif a <= evalf(2*e*sin((1/2)*B)) then eq2(a) 
  else 0 
  end if 
end proc:

@Heeka 

(1) Y_vals produces a list of y(x), for x=0, 0.1, 0.2, ..., 1.0.

(2) Comparison exact and numeric solution:

restart;
DE := (x+b*y(x))*diff(y(x),x) + y(x) = 0;
BC := y(1)=1:
b := 1/10:
# exact solution:
sol1 := dsolve( {DE,BC}, y(x) ): Y1 := unapply( subs( sol1, y(x) ), x );
# numeric solution:
sol := dsolve( {DE,BC}, y(x), numeric, output=listprocedure ):
Y := subs( sol, y(x) ):
# comparison exact and numeric sulution
plot( Y1-Y, 0..1 );