#------- Utilities

`diff/piecewise` := proc ()
  description  `cheap diff on piecewise expression`, 
               `Maple 7's library function considers`,
               `discontinuities, but is _much_ slower.`;
  local g, s, n, k, t;
  g := [args[1..nargs-1]];
  t := args[nargs];  n := nops(g);
  s := seq(op([g[2*k-1],diff(g[2*k],t)]),k=1..floor(n/2));
  if (type(n,odd)) then s := s,diff(g[n],t); fi;
  RETURN (piecewise(s));
end:

signum_rightof := proc (f, t::numeric)
  description `return signum of f(t+)`;
  global eps;
  local r, s, x;
  if (type (t, rational)) then
    s := simplify (series (f(x), x=t));
    r := op(1,s);
    RETURN (signum(r));
  fi;
  RETURN (signum(evalf (f(t+eps))));
end:

fsolve_smallest := proc (f, r::range)
  description `return smallest zero in [r]. NULL if none.`;
  global eps;
  local tfrom, tto, s, sprev, t;
  tfrom := op(1,r); tto := op(2,r);
  sprev := NULL;
  while (tfrom <= tto) do
    s := fsolve (f(t), t, t=tfrom..tto);
    if (whattype(s) <> float) then RETURN(sprev); fi;
    if ((s < tfrom) or (tto < s)) then RETURN(sprev); fi;
    tto := s-eps;
    sprev := s;
  od;
  RETURN(sprev);
end:

#------- get next interval

case[1,0] := 1:  case[0,1] := 2:  case[1,1] := 3:
# i[2] = 0:  i[1] = 0:

next_case := proc (t0::numeric, j::array)
  description `returns case that follows point t0.`;
  global eps, case;
  RETURN (case[signum_rightof(j[1],t0), signum_rightof(j[2],t0)]);
end:

next_interval := proc (t0::numeric, i0::array, VAR_case)
  description `returns currents for next interval.`,
              `In VAR_case next case.`;
  local t, j, s, case;
  for s from nof_cases by -1 to 1 do
    j := array(1..2);
    j[1] := unapply (i[s][1](t,t0,i0[1],i0[2]), t);
    j[2] := unapply (i[s][2](t,t0,i0[1],i0[2]), t);
    case := next_case (t0, j);
    if (s = case) then VAR_case := s; RETURN (j); fi;
  od;
  ERROR (`Continuation failed`);
end:

crit_smallest := proc (f, r::range(numeric))
  description `smallest zero of f in ]r] (without left border).`,
              `right(r) if (f == 0) or (f(t) <> 0) in r.`;
  global eps;
  local s, tfrom, tto;
  tfrom := op(1,r); tto := op(2,r);
  if (f = 0) then RETURN (tto); fi;
  s := fsolve_smallest (f, tfrom+eps..tto);
  if (whattype(s) <> float) then RETURN (tto); fi;
  RETURN (s);
end:

next_crit := proc (j::array, r::range(numeric), 
                   case::integer, VAR_i)
  description `returns next critical point.`,
              `In VAR_i current i[] at this point.`;
  local t, ug, ud, crit, tcrit, i;
  ug := unapply (evalf(R3*(j[1](t) + j[2](t)) +
        L3*(diff(j[1](t),t)+diff(j[2](t),t))), t);
  ud[1] := unapply (u[1](t) - ug(t), t);
  ud[2] := unapply (u[2](t) - ug(t), t);
  if (case = 1) then
    crit[1] := crit_smallest (j[1], r);
    crit[4] := crit_smallest (ud[2], r);
    tcrit := min(crit[1], crit[4]);
  elif (case = 2) then
    crit[2] := crit_smallest (j[2], r);
    crit[3] := crit_smallest (ud[1], r);
    tcrit := min(crit[2], crit[3]);
  elif (case = 3) then
    crit[1] := crit_smallest (j[1], r);
    crit[2] := crit_smallest (j[2], r);
    crit[3] := crit_smallest (ud[1], r);
    crit[4] := crit_smallest (ud[2], r);
    tcrit := min(crit[1], crit[2], crit[3], crit[4]);
  fi;
  i := array (1..2, [j[1](tcrit), j[2](tcrit)]);
  if (crit[1] = tcrit) then i[1] := 0; fi;
  if (crit[2] = tcrit) then i[2] := 0; fi;
  VAR_i := i;
  RETURN (tcrit);
end:

#------- main program

rectify_solve := proc (i0::array(numeric), r::range(numeric))
  description `returns i[1..2] as piecewise functions`;
  global nof_states, eps;
  local tfrom, tto, tcur, scur, res, i, j, k, ires, n, nof, d, s, t;
  tfrom := op(1,r); tto := op(2,r);
  res := array (1..nof_states);
  for k from 1 to nof_states do res[k] := NULL; od;
  tcur := tfrom; i := i0;
  while (tcur < tto) do
    unassign ('scur');
    j := next_interval (tcur, i, scur);
    for k from 1 to nof_states do 
      res[k] := res[k], [tcur, j[k], scur]; 
    od;
    unassign ('i');
    tcur := next_crit (j, tcur..tto, scur, i);
  od;
  ires := array (1..nof_states);
  for k from 1 to nof_states do 
    nof := nops([res[k]]);
    d := res[k],[tto+eps];
    s := seq(op([d[n][1] <= t and t < d[n+1][1],d[n][2](t)]),n=1..nof);
    ires[k] := unapply(piecewise(s),t);
  od;
  RETURN (ires);
end: