I did a bit of experimenting, and it seems that the following example code does what I need:

a[1]=0;
a[2]=1;
a[3]=3;
b[1]=1;
b[2]=3;
b[3]=6;
f[1][x_]:=10;
f[2][x_]:=20;
f[3][x_]:=30;
nintervals=3;
Fun[x_]:=Module[{n,v},
If[x<a[1] || x>b[nintervals],v=Indeterminate,For[n=1,n<=nintervals,n=n+1, If[x>=a[n]&&x<=b[n],v=f[n][x];Break[]]]];
v];

The trick is to use the Break[] command, to interrupt the loop. I did not realise it was possible, as the help files provide only the two variants of If[]:

If[condition,t,f]

and

If[condition,t,f,u]

but don't provide

If[condition,t]

which is needed here.

I have tried the version with Piecewise + Table, but in such a case a get a strange error. In my program, functions f[n][] are calculated by Interpolating Function objects obj[1] and obj[2] obtained by NDSolve, while solving some ODEs. My code is:

obj[1]=Import["obj1.m"];
obj[2]=Import["obj2.m"];
nmax=2;
For[n=1,n<=nmax,n=n+1,heads[n]=Head[(Subscript[psix, 1][th]/.obj[n])]];
For[n=1,n<=nmax,n=n+1,a[n]=SetPrecision[First[First[heads[n]["Domain"]]],Infinity]];
For[n=1,n<=nmax,n=n+1,b[n]=SetPrecision[Last[First[heads[n]["Domain"]]],Infinity]];
bounds=N[Table[{a[n],b[n]},{n,1,nmax,1}],25]
{{881.3511598761730369258651,52027.18417423114383137352},{52027.18417423114383137352,1.093770863505275159944475*10^7}}
For[n=1,n<=nmax,n=n+1,f[n][th_]:=Evaluate[(Subscript[psix, 1][th]/.obj[n])]];
Fun1[th_]:=Piecewise[Table[{f[n][th],th>a[n]&&th<=b[n]},{n,1,nmax}],"Undefined"];
Fun2[th_]:=Module[{n,v},
If[th<a[1] || th>b[nmax],v=Indeterminate,For[n=1,n<=nmax,n=n+1, If[th>=a[n]&&th<=b[n],v=f[n][th];Break[]]]];
v];
Fun1[10001.176032248239`]
InterpolatingFunction::dmval: Input value {10001.2} lies outside the range of data in the interpolating function. Extrapolation will be used.
1.11039
f[1][10001.176032248239`]
1.11039
Fun2[10001.176032248239`]
1.11039

As can be seen, functions Fun1[th], f[1][th] and Fun2[th] all return the same (seemingly correct) value for a certain argument, but in the case of Fun1[th], which uses Piecewise[], there is an incomprehensible error message. It is incomprehensible, because the argument value lies within the domain covered by obj[1]. How to understand what happens here? Unfortunately I cannot attach the obj1.m and obj2.m files, because they have ca. 40MB each.

Lesław

You have some unnecessary and overly complicated definitions. I'd suggest something along the following lines (I don't know your whole context, so there are probably further improvements, but hopefully this will get you started):

First, encode your rules as data:

intervalFunctionMap = 
  <|Interval[{0, 1}] -> f[1], 
    Interval[{1, 3}] -> f[2], 
    Interval[{3, 6}] -> f[3]|>

Now just define your function as a selection:

Fun[x_] := 
  First[KeySelect[intervalFunctionMap, IntervalMemberQ[#, x] &], Indeterminate &][x]

Demo:

Fun /@ Range[-1, 7, .5]
(* {Indeterminate, Indeterminate, f[1][0.], f[1][0.5], f[1][1.], 
    f[2][1.5], f[2][2.], f[2][2.5], f[2][3.], f[3][3.5], f[3][4.], 
    f[3][4.5], f[3][5.], f[3][5.5], 
    f[3][6.], Indeterminate, Indeterminate} *)

Piecewise is a good choice for this. Note that Piecewise evaluates the conditions in turn, until one of them evaluates to True. If the intervals are contiguous and listed in order, then only one boundary need be specified for each. A new interval could be added to the chain as long as it is added at the correct place in the order.

Like this:

f[1] = 1; f[2] = 2; f[3] = 3;

f[x_] := Piecewise[{
   {f[1], x < 1},
   {f[2], x < 2},
   {f[3], x < 3},
   {5, True}
   }]

Plot[f[x], {x, 0, 5}]

The f[n] can of course be functions of x rather than constants.