# HG changeset patch
# User Carlo de Falco <carlo.defalco@polimi.it>
# Date 1460176613 -7200
# Node ID 07d30e6fcfdee01d9146d580d3500b9c14351863
# Parent  3933d5073baabc4ae8b857fce0d7d878eae65fc1
Add demos to ode23 and ode45

* scripts/ode/ode{23,45}.m : Add new demo demonstrating
  order of convergence.

diff -r 3933d5073baa -r 07d30e6fcfde scripts/ode/ode23.m
--- a/scripts/ode/ode23.m	Fri Apr 08 16:51:05 2016 -0400
+++ b/scripts/ode/ode23.m	Sat Apr 09 06:36:53 2016 +0200
@@ -590,3 +590,28 @@
 %!  opt = struct ("foo", 1);
 %!  [t, y] = ode23 (@fpol, [0 2], [2 0], opt);
 
+%!demo
+%!
+%! # Demonstrate convergence order for ode23
+%! tol = 1e-5 ./ 10.^[0:8];
+%! for ii = 1 : numel (tol)
+%!   opt = odeset ("RelTol", tol(ii), "AbsTol", realmin);
+%!   [t, y] = ode23 (@(t, y) -y, [0, 1], 1, opt);
+%!   h(ii) = 1 / (numel (t) - 1);
+%!   err(ii) = norm (y .* exp (t) - 1, inf);
+%! endfor
+%!
+%! # Estimate order numerically
+%! p = diff (log (err)) ./ diff (log (h))
+%!
+%! %! # Estimate order visually
+%! figure;
+%! loglog (h, tol, '-ob', h, err,
+%!         h, (h/h(end)) .^ 2 .* tol(end), "k--",
+%!         h, (h/h(end)) .^ 3 .* tol(end), "k-")
+%! axis tight
+%! xlabel ("h");
+%! ylabel ("err(h)");
+%! title ("Convergence plot for ode23");
+%! legend ("imposed tolerance", "ode23 (relative) error",
+%!         "order 2", "order 3");
diff -r 3933d5073baa -r 07d30e6fcfde scripts/ode/ode45.m
--- a/scripts/ode/ode45.m	Fri Apr 08 16:51:05 2016 -0400
+++ b/scripts/ode/ode45.m	Sat Apr 09 06:36:53 2016 +0200
@@ -587,3 +587,29 @@
 %! opt = struct ("foo", 1);
 %! [t, y] = ode45 (@fpol, [0 2], [2 0], opt);
 
+
+%!demo
+%!
+%! # Demonstrate convergence order for ode45
+%! tol = 1e-5 ./ 10.^[0:8];
+%! for ii = 1 : numel (tol)
+%!   opt = odeset ("RelTol", tol(ii), "AbsTol", realmin);
+%!   [t, y] = ode45 (@(t, y) -y, [0, 1], 1, opt);
+%!   h(ii) = 1 / (numel (t) - 1);
+%!   err(ii) = norm (y .* exp (t) - 1, inf);
+%! endfor
+%!
+%! # Estimate order numerically
+%! p = diff (log (err)) ./ diff (log (h))
+%!
+%! %! # Estimate order visually
+%! figure;
+%! loglog (h, tol, '-ob', h, err,
+%!         h, (h/h(end)) .^ 4 .* tol(end), "k--",
+%!         h, (h/h(end)) .^ 5 .* tol(end), "k-")
+%! axis tight
+%! xlabel ("h");
+%! ylabel ("err(h)");
+%! title ("Convergence plot for ode45");
+%! legend ("imposed tolerance", "ode45 (relative) error",
+%!         "order 4", "order 5");