Graphical Solution

3D plot (note that Maple is able to handle epsilon=0 in a reasonable way

> plot3d( rhs(SOLN), t=0..1, epsilon=0..0.1 );

>

2D plots for selected values of epsilon. Note the rapid initial growth that occurs as [Maple Math] -> 0; the solution is nearly linear closer to 1.

> plot( [ seq( eval( rhs(SOLN), epsilon=e ), e=[0.001,0.01,0.1,1] )], t=0..1 );