There is a bug in example 4, which cause it to not follow the boundary values,
the correct res. function should be:
res = [ ya(1)
yb(1)
(ya(2)-yb(2)) ];
which is not dependent on T, instead of
res = [ya(1) - yb(1)
ya(2) - yb(2)
T*(-1/3)*(ya(1) - 0.7 + 0.8*ya(2)) - 1];
Comment only