function odeshow;

   ProblemsList = {'None', [];
                   'Euler for stiff problem', @StiffExample;
                   'Example with singular point', @SingularPointExample};

   % Default values:
   ExplicitEulerRadioInitValue = 1;
   ExplicitEulerRightRadioInitValue = 1;
   ExplicitEulerLeftRadioInitValue = 0;
   ExplicitEulerStepEditInitValue = num2str(.1);
   ExplicitEulerLocalErrorRadioInitValue = 0;
   ImplicitEulerRadioInitValue = 0;
   ImplicitEulerRightRadioInitValue = 0;
   ImplicitEulerLeftRadioInitValue = 0;
   ImplicitEulerStepEditInitValue = num2str(.1);
   DirectionsRadioInitValue = 0;
   tDirectionsPointsNumberEditInitValue = num2str(10);
   yDirectionsPointsNumberEditInitValue = num2str(10);
   ode45RadioInitValue = 0;
   ode15sRadioInitValue = 0;
   ode23sRadioInitValue = 0;
   fEditInitValue = '-10.*t.*y';
   aEditInitValue = num2str(-1);
   bEditInitValue = num2str(-.4);
   y0EditInitValue = num2str(.05);
   ProblemsPopUpInitValue = 1;
   YLimitsRadioInitValue = 0;
   YLimitsMinEditInitValue = 0;
   YLimitsMaxEditInitValue = 3.5;

   % Figure:
   Fig = openfig('odeshow.fig');
   handles = guihandles(Fig);
   
   % Its uicontrols:
   Axes = handles.Axes;
   CloseButton = handles.CloseButton;
   GoButton = handles.GoButton;
   ImplicitEulerRadio = handles.ImplicitEulerRadio;
   ImplicitEulerRightRadio = handles.ImplicitEulerRightRadio;
   ImplicitEulerLeftRadio = handles.ImplicitEulerLeftRadio;
   ImplicitEulerStepEdit = handles.ImplicitEulerStepEdit;
   ExplicitEulerRadio = handles.ExplicitEulerRadio;
   ExplicitEulerRightRadio = handles.ExplicitEulerRightRadio;
   ExplicitEulerLeftRadio = handles.ExplicitEulerLeftRadio;
   ExplicitEulerStepEdit = handles.ExplicitEulerStepEdit;
   ExplicitEulerLocalErrorRadio = handles.ExplicitEulerLocalErrorRadio;
   DirectionsRadio = handles.DirectionsRadio;
   tDirectionsPointsNumberEdit = handles.tDirectionsPointsNumberEdit;
   yDirectionsPointsNumberEdit = handles.yDirectionsPointsNumberEdit;
   ode45Radio = handles.ode45Radio;
   ode15sRadio = handles.ode15sRadio;
   ode23sRadio = handles.ode23sRadio;
   ode45Text = handles.ode45Text;
   ode15sText = handles.ode15sText;
   ode23sText = handles.ode23sText;
   fEdit = handles.fEdit;
   aEdit = handles.aEdit;
   bEdit = handles.bEdit;
   y0Edit = handles.y0Edit;
   ProblemsPopUp = handles.ProblemsPopUp;
   YLimitsRadio = handles.YLimitsRadio;
   YLimitsMinEdit = handles.YLimitsMinEdit;
   YLimitsMaxEdit = handles.YLimitsMaxEdit;

   % Init values:
   set(Fig, 'DoubleBuffer', 'On');
   set(Axes, 'XGrid', 'On', 'YGrid','On', 'NextPlot', 'Add');
   set(CloseButton, 'CallBack', 'close');
   set(GoButton, 'CallBack', @GoButtonCallBack);
   set(ProblemsPopUp, ...
     'Value', ProblemsPopUpInitValue, ...
     'String', ProblemsList(:, 1), ...
     'CallBack', @ProblemsPopUpCallBack);
   SetInitValues;
   
   % Plots:
   lst = {0, 0, 'Parent', Axes, 'Visible', 'off', 'Color'};
   ExplicitEulerPlot = line(lst{:}, 'r', 'Marker', 'o');
   ExplicitEulerLeftPlot = line(lst{:}, 'b', 'LineStyle', '--');
   ExplicitEulerRightPlot = line(lst{:}, 'b', 'LineStyle', '--');
   ExplicitEulerLocalErrorPlot = line(lst{:}, 'm');
   ImplicitEulerPlot = line(lst{:}, 'k', 'Marker', 'o');
   ImplicitEulerLeftPlot = line(lst{:}, 'b', 'LineStyle', '--');
   ImplicitEulerRightPlot = line(lst{:}, 'b', 'LineStyle', '--');
   ode45Plot = line(lst{:}, 'm', 'Marker', 's');
   ode15sPlot = line(lst{:}, 'g', 'Marker', 'd');
   ode23sPlot = line(lst{:}, 'c', 'Marker', 'p');
   DirectionsPlot = line(lst{:}, 'k');

   RecomputeAndRedraw;

   function SetInitValues
       set(ImplicitEulerRadio, ...
         'Value', ImplicitEulerRadioInitValue, ...
         'CallBack', @ImplicitEulerRadioCallBack);
       set(ImplicitEulerRightRadio, ...
         'Value', ImplicitEulerRightRadioInitValue, ...
         'CallBack', @ImplicitEulerRightRadioCallBack);
       set(ImplicitEulerLeftRadio, ...
         'Value', ImplicitEulerLeftRadioInitValue, ...
         'CallBack', @ImplicitEulerLeftRadioCallBack);
       set(ImplicitEulerStepEdit, ...
         'String', ImplicitEulerStepEditInitValue, ...
         'CallBack', @ImplicitEulerStepEditCallBack);
       set(ExplicitEulerRadio, ...
         'Value', ExplicitEulerRadioInitValue, ...
         'CallBack', @ExplicitEulerRadioCallBack);
       set(ExplicitEulerRightRadio, ...
         'Value', ExplicitEulerRightRadioInitValue, ...
         'CallBack', @ExplicitEulerRightRadioCallBack);
       set(ExplicitEulerLeftRadio, ...
         'Value', ExplicitEulerLeftRadioInitValue, ...
         'CallBack', @ExplicitEulerLeftRadioCallBack);
       set(ExplicitEulerStepEdit, ...
         'String', ExplicitEulerStepEditInitValue, ...
         'CallBack', @ExplicitEulerStepEditCallBack);
       set(ExplicitEulerLocalErrorRadio, ...
         'Value', ExplicitEulerLocalErrorRadioInitValue, ...
         'CallBack', @ExplicitEulerLocalErrorRadioCallBack);
       set(DirectionsRadio, ...
         'Value', DirectionsRadioInitValue, ...
         'CallBack', @DirectionsRadioCallBack);
       set(tDirectionsPointsNumberEdit, ...
         'String', tDirectionsPointsNumberEditInitValue, ...
         'CallBack', @tDirectionsPointsNumberEditCallBack);
       set(yDirectionsPointsNumberEdit, ...
         'String', yDirectionsPointsNumberEditInitValue, ...
         'CallBack', @yDirectionsPointsNumberEditCallBack);
       set(ode45Radio, ...
         'Value', ode45RadioInitValue, ...
         'CallBack', @ode45RadioCallBack);
       set(ode23sRadio, ...
         'Value', ode23sRadioInitValue, ...
         'CallBack', @ode23sRadioCallBack);
       set(ode15sRadio, ...
         'Value', ode15sRadioInitValue, ...
         'CallBack', @ode15sRadioCallBack);
       set(fEdit, ...
         'String', fEditInitValue, ...
         'CallBack', @fEditCallBack);
       set(aEdit, ...
         'String', aEditInitValue, ...
         'CallBack', @aEditCallBack);
       set(bEdit, ...
         'String', bEditInitValue, ...
         'CallBack', @bEditCallBack);
       set(y0Edit, ...
         'String', y0EditInitValue, ...
         'CallBack', @y0EditCallBack);
       set(YLimitsRadio, ...
         'Value', YLimitsRadioInitValue, ...
         'CallBack', @YLimitsRadioCallBack);
       set(YLimitsMinEdit, ...
         'String', YLimitsMinEditInitValue, ...
         'CallBack', @YLimitsMinEditCallBack);
       set(YLimitsMaxEdit, ...
         'String', YLimitsMaxEditInitValue, ...
         'CallBack', @YLimitsMaxEditCallBack);
   end;

   function ExplicitEulerRadioCallBack(varargin)
     if get(ExplicitEulerRadio, 'Value') == 1
       set(ExplicitEulerPlot, 'Visible', 'on');
       set([ExplicitEulerRightRadio, ExplicitEulerLeftRadio, ...
            ExplicitEulerLocalErrorRadio], 'Enable', 'on');
       ExplicitEulerRightRadioCallBack;
       ExplicitEulerLeftRadioCallBack;
       ExplicitEulerLocalErrorRadioCallBack;
     else
       set([ExplicitEulerPlot, ExplicitEulerRightPlot, ...
            ExplicitEulerLeftPlot, ExplicitEulerLocalErrorPlot], 'Visible', 'off');
       set([ExplicitEulerRightRadio, ExplicitEulerLeftRadio, ...
            ExplicitEulerLocalErrorRadio], 'Enable', 'off');
     end;
   end;

   function ExplicitEulerRightRadioCallBack(varargin)
     if get(ExplicitEulerRadio, 'Value') == 1 && get(ExplicitEulerRightRadio, 'Value') == 1
       set(ExplicitEulerRightPlot, 'Visible', 'on');
     else
       set(ExplicitEulerRightPlot, 'Visible', 'off');
     end;
   end;

   function ExplicitEulerLeftRadioCallBack(varargin)
     if get(ExplicitEulerRadio, 'Value') == 1 && get(ExplicitEulerLeftRadio, 'Value') == 1
       set(ExplicitEulerLeftPlot, 'Visible', 'on');
     else
       set(ExplicitEulerLeftPlot, 'Visible', 'off');
     end;
   end;

   function ExplicitEulerLocalErrorRadioCallBack(varargin)
     if get(ExplicitEulerLocalErrorRadio, 'Value') == 1 && get(ExplicitEulerRadio, 'Value') == 1
       set(ExplicitEulerLocalErrorPlot, 'Visible', 'on');
     else
       set(ExplicitEulerLocalErrorPlot, 'Visible', 'off');
     end;
   end;

   function GoButtonCallBack(varargin)
      RecomputeAndRedraw;
   end;

   function ExplicitEulerStepEditCallBack(varargin)
   end;

   function ImplicitEulerRadioCallBack(varargin)
     if get(ImplicitEulerRadio, 'Value') == 1
       set(ImplicitEulerPlot, 'Visible', 'on');
       set([ImplicitEulerRightRadio, ImplicitEulerLeftRadio], 'Enable', 'on');
       ImplicitEulerRightRadioCallBack;
       ImplicitEulerLeftRadioCallBack;
     else
       set([ImplicitEulerPlot, ImplicitEulerRightPlot, ImplicitEulerLeftPlot], 'Visible', 'off');
       set([ImplicitEulerRightRadio, ImplicitEulerLeftRadio], 'Enable', 'off');
     end;
   end;

   function ImplicitEulerRightRadioCallBack(varargin)
     if get(ImplicitEulerRadio, 'Value') == 1 && get(ImplicitEulerRightRadio, 'Value') == 1
       set(ImplicitEulerRightPlot, 'Visible', 'on');
     else
       set(ImplicitEulerRightPlot, 'Visible', 'off');
     end;
   end;

   function ImplicitEulerLeftRadioCallBack(varargin)
     if get(ImplicitEulerRadio, 'Value') == 1 && get(ImplicitEulerLeftRadio, 'Value') == 1
       set(ImplicitEulerLeftPlot, 'Visible', 'on');
     else
       set(ImplicitEulerLeftPlot, 'Visible', 'off');
     end;
   end;

   function ImplicitEulerStepEditCallBack(varargin)
   end;

   function DirectionsRadioCallBack(varargin)
     if get(DirectionsRadio, 'Value') == 1
       set(DirectionsPlot, 'Visible', 'On');
     else
       set(DirectionsPlot, 'Visible', 'Off');
     end;
   end;

   function tDirectionsPointsNumberEditCallBack(varargin)
      RecomputeAndRedrawDirections;
   end;

   function yDirectionsPointsNumberEditCallBack(varargin)
      RecomputeAndRedrawDirections;
   end;

   function ode45RadioCallBack(varargin)
     if get(ode45Radio, 'Value') == 1
       set(ode45Plot, 'Visible', 'on');
     else
       set(ode45Plot, 'Visible', 'off');
     end;
   end;

   function ode15sRadioCallBack(varargin)
     if get(ode15sRadio, 'Value') == 1
       set(ode15sPlot, 'Visible', 'on');
     else
       set(ode15sPlot, 'Visible', 'off');
     end;
   end;

   function ode23sRadioCallBack(varargin)
     if get(ode23sRadio, 'Value') == 1
       set(ode23sPlot, 'Visible', 'on');
     else
       set(ode23sPlot, 'Visible', 'off');
     end;
   end;

   function fEditCallBack(varargin)
   end;

   function aEditCallBack(varargin)
   end;

   function bEditCallBack(varargin)
   end;

   function y0EditCallBack(varargin)
   end;

   function ProblemsPopUpCallBack(varargin)
     no = get(ProblemsPopUp, 'Value');
     if ~isempty(ProblemsList{no, 2})
       feval(ProblemsList{no, 2});
       SetInitValues;
     end
   end;

   function YLimitsRadioCallBack(varargin)
      if get(YLimitsRadio, 'Value') == 1
        set([YLimitsMaxEdit, YLimitsMinEdit], 'Enable', 'on');
        SetAxesYLimits;
      else
        DeleteDirections;
        set(Axes, 'YLimMode', 'Auto');
        set([YLimitsMaxEdit, YLimitsMinEdit], 'Enable', 'off');
        RecomputeAndRedrawDirections;
      end;
   end;

   function YLimitsMaxEditCallBack(varargin)
      SetAxesYLimits;
   end;

   function YLimitsMinEditCallBack(varargin)
      SetAxesYLimits;
   end;
   
   function SetAxesYLimits
        y_min = str2num(get(YLimitsMinEdit, 'String'));
        y_max = str2num(get(YLimitsMaxEdit, 'String'));
        set(Axes, 'YLim', [y_min, y_max]);
        RecomputeAndRedrawDirections;
   end;

   function RecomputeAndRedrawDirections
     f = inline([get(fEdit, 'String'), ' + 0*y + 0*t']);
     t = get(Axes, 'XLim');
     y = get(Axes, 'YLim');
     tn = str2num(get(tDirectionsPointsNumberEdit, 'String')) + 1;
     yn = str2num(get(yDirectionsPointsNumberEdit, 'String')) + 1;
     tt = linspace(t(1), t(2), tn);
     yy = linspace(y(1), y(2), yn);
     [T, Y] = meshgrid(tt, yy);
     F = feval(f, T, Y);
     alpha = 20; %pixels
     DT = ones(size(F));
     DY = F;
     set(Axes, 'Units', 'Pixels');
     P = get(Axes, 'Position');
     Len = sqrt((DT .* P(3) ./ (t(2) - t(1))).^2 + ...
                (DY .* P(4) ./ (y(2) - y(1))).^2);
     [T, Y] = directions(T, Y, alpha * DT ./ Len, alpha * DY ./Len);
     set(DirectionsPlot, 'XData', T, 'YData', Y);
   end;

   function DeleteDirections
     set(DirectionsPlot, 'XData', NaN, 'YData', NaN);
   end;
   
   function RecomputeAndRedraw;
     Recompute;
     Redraw;
     RecomputeAndRedrawDirections;
   end;

   function Recompute;
   
     a = str2num(get(aEdit, 'String'));
     b = str2num(get(bEdit, 'String'));
     if a > b
        msgbox('a should be less than b', 'Error', 'error', 'modal');
        return;
     end  
     y0 = str2num(get(y0Edit, 'String'));
     f = inline([get(fEdit, 'String'), ' + 0*y + 0*t']);
     h_expl = str2num(get(ExplicitEulerStepEdit, 'String'));
     h_impl = str2num(get(ImplicitEulerStepEdit, 'String'));

     % Далее в решателях иногда можно использовать f, а иногда только @fn
     % Описание fn см. ниже
     [t_expl, y_expl] = euler(f, [a, b], y0, h_expl); 
     [t_impl, y_impl] = eulerimplicit(f, [a, b], y0, h_impl);
     [t_ode45, y_ode45] = ode45(f, [a, b], y0);
     [t_ode15s, y_ode15s] = ode15s(f, [a, b], y0);
     [t_ode23s, y_ode23s] = ode23s(f, [a, b], y0);

     set(ode45Text, 'String', num2str(length(t_ode45)));
     set(ode15sText, 'String', num2str(length(t_ode15s)));
     set(ode23sText, 'String', num2str(length(t_ode23s)));

     set(ExplicitEulerPlot, 'XData', t_expl, 'YData', y_expl);
     set(ImplicitEulerPlot, 'XData', t_impl, 'YData', y_impl);
     set(ode45Plot, 'XData', t_ode45, 'YData', y_ode45);
     set(ode15sPlot, 'XData', t_ode15s, 'YData', y_ode15s);
     set(ode23sPlot, 'XData', t_ode23s, 'YData', y_ode23s);

     XLim = [a, b];
     set(Axes, 'XLim', XLim);

     t_expl_right = [];
     y_expl_right = [];
     t_expl_left = [];
     y_expl_left = [];
     t_loc_err = [];
     y_loc_err = [];
     n = length(t_expl);
     for k = 1:n;
       if k ~= n && t_expl(k) < XLim(2) - eps % Иначе зависание
         sol = ode23s(@fn, [t_expl(k), XLim(2)], y_expl(k));
         t_expl_right = [t_expl_right, sol.x, NaN];
         y_expl_right = [y_expl_right, sol.y, NaN];
         if t_expl(k + 1) <= sol.x(end)
           t_loc_err = [t_loc_err, t_expl(k + 1), t_expl(k + 1), NaN];
           y_loc_err = [y_loc_err, y_expl(k + 1), deval(sol, t_expl(k + 1)), NaN];
         end;
       end;
       if t_expl(k) > XLim(1) + eps % Иначе зависание
         [t, y] = ode23s(f, [t_expl(k), XLim(1)], y_expl(k));
         t_expl_left = [t_expl_left, t', NaN];
         y_expl_left = [y_expl_left, y', NaN];
       end;
     end;
     set(ExplicitEulerRightPlot, 'XData', t_expl_right, 'YData', y_expl_right);
     set(ExplicitEulerLeftPlot, 'XData', t_expl_left, 'YData', y_expl_left);
     set(ExplicitEulerLocalErrorPlot, 'XData', t_loc_err, 'YData', y_loc_err);

     t_impl_right = [];
     y_impl_right = [];
     t_impl_left = [];
     y_impl_left = [];
     n = length(t_impl);
     for k = 1:n;
       if k ~= n && t_impl(k) < XLim(2) - eps % Иначе зависание
         [t, y] = ode23s(f, [t_impl(k), XLim(2)], y_impl(k));
         t_impl_right = [t_impl_right, t', NaN];
         y_impl_right = [y_impl_right, y', NaN];
       end;
       if t_impl(k) > XLim(1) + eps % Иначе зависание
         [t, y] = ode23s(f, [t_impl(k), XLim(1)], y_impl(k));
         t_impl_left = [t_impl_left, t', NaN];
         y_impl_left = [y_impl_left, y', NaN];
       end;
     end;
     set(ImplicitEulerRightPlot, 'XData', t_impl_right, 'YData', y_impl_right);
     set(ImplicitEulerLeftPlot, 'XData', t_impl_left, 'YData', y_impl_left);

     function v = fn(t, y, varargin)
       v = feval(f, t, y, varargin{:});
     end;
     
   end;

   function Redraw
     ode45RadioCallBack;
     ode15sRadioCallBack;
     ode23sRadioCallBack;
     ExplicitEulerRadioCallBack;
     ExplicitEulerLeftRadioCallBack;
     ExplicitEulerRightRadioCallBack;
     ImplicitEulerRadioCallBack;
     ImplicitEulerLeftRadioCallBack;
     ImplicitEulerRightRadioCallBack;
     YLimitsRadioCallBack;
     DirectionsRadioCallBack;
   end;

   function [X, Y] = directions(x, y, u, v)
     n = length(x(:));
     X = [x(:), x(:) + u(:), NaN*zeros(n, 1)]';
     X = X(:);
     Y = [y(:), y(:) + v(:), NaN*zeros(n, 1)]';
     Y = Y(:);
   end

   function StiffExample
      ExplicitEulerRadioInitValue = 1;
      ExplicitEulerRightRadioInitValue = 1;
      ExplicitEulerLeftRadioInitValue = 1;
      ExplicitEulerStepEditInitValue = num2str(.1);
      ExplicitEulerLocalErrorRadioInitValue = 0;
      ImplicitEulerRadioInitValue = 0;
      ImplicitEulerRightRadioInitValue = 0;
      ImplicitEulerLeftRadioInitValue = 0;
      ImplicitEulerStepEditInitValue = num2str(.1);
      DirectionsRadioInitValue = 0;
      tDirectionsPointsNumberEditInitValue = num2str(10);
      yDirectionsPointsNumberEditInitValue = num2str(10);
      ode45RadioInitValue = 0;
      ode15sRadioInitValue = 0;
      ode23sRadioInitValue = 0;
      fEditInitValue = '-21*(y-sin(t))+cos(t)';
      aEditInitValue = num2str(0);
      bEditInitValue = num2str(1);
      y0EditInitValue = num2str(1);
      YLimitsRadioInitValue = 1;
      YLimitsMinEditInitValue = -2;
      YLimitsMaxEditInitValue = 4;
   end;

   function SingularPointExample
      ExplicitEulerRadioInitValue = 1;
      ExplicitEulerRightRadioInitValue = 0;
      ExplicitEulerLeftRadioInitValue = 0;
      ExplicitEulerStepEditInitValue = num2str(.1);
      ExplicitEulerLocalErrorRadioInitValue = 0;
      ImplicitEulerRadioInitValue = 1;
      ImplicitEulerRightRadioInitValue = 0;
      ImplicitEulerLeftRadioInitValue = 0;
      ImplicitEulerStepEditInitValue = num2str(.1);
      DirectionsRadioInitValue = 1;
      tDirectionsPointsNumberEditInitValue = num2str(10);
      yDirectionsPointsNumberEditInitValue = num2str(10);
      ode45RadioInitValue = 1;
      ode15sRadioInitValue = 1;
      ode23sRadioInitValue = 1;
      fEditInitValue = 'y./(3*t)';
      aEditInitValue = num2str(-1);
      bEditInitValue = num2str(1);
      y0EditInitValue = num2str(-1);
      YLimitsRadioInitValue = 1;
      YLimitsMinEditInitValue = -1;
      YLimitsMaxEditInitValue = 1;
   end;

end