readlib(readdata);
ucalc := readdata(`tmp.u`,2):
upcalc := readdata(`tmp.up`,2):
a := 2; 
b := 1;		# choose various cases for testing
exa := x -> (a+x^(1/2))*(1-x)^(3/2)+b*x;	# exact answer, for checking
ep := unapply(diff(exa(x),x),x);

pld := proc(data, func)  # plot function at same abscissa as given data set
   local ret,i;
   ret := [];
   for i in data do;
      ret := [op(ret), [i[1], func(i[1])]];
   od;
   ret;
end;

# compare abscissas and ordinates of given 
# data sets, and return the difference in ordinates for those
# entries with matching abscissas  (dat2 -dat1)
compa := proc(dat1, dat2)
   local i, ret, d1, d2;
   ret := [];
   for i from 1 to nops(dat1) do;
      d1 := dat1[i]; d2 := dat2[i];
      if (d1[1] = d2[1]) then;
	 ret := [op(ret), [d1[1], d2[2]- d1[2]]];
      fi;
   od;
   ret;
end;
#	as compa, but relative comparison: (dat2-dat1)/dat2   
compr := proc(dat1, dat2)
   local i, ret, d1, d2;
   ret := [];
   for i from 1 to nops(dat1) do;
      d1 := dat1[i]; d2 := dat2[i];
      if (d1[1] = d2[1]) then;
	 ret := [op(ret), [d1[1], (d2[2]- d1[2])/d2[2]]];
      fi;
   od;
   ret;
end;


pl := proc(f,a,b)
   local i, t;
   t := (b-a)/25;
   [seq(evalf([a+t*i, f(a+t*i)]), i=0..25)];
end;