From ba154db78b02163be343b235438e3add73cba526 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Fri, 14 Sep 2012 23:51:30 -0400
Subject: [PATCH] Rewrite the forward_euler() function to return a vector of
 coefficients rather than assuming that we will really have a function to
 evaluate.

---
 forward_euler.m | 44 +++++++++++++++++++++-----------------------
 1 file changed, 21 insertions(+), 23 deletions(-)

diff --git a/forward_euler.m b/forward_euler.m
index c6ef217..ec678a1 100644
--- a/forward_euler.m
+++ b/forward_euler.m
@@ -1,26 +1,20 @@
-function df = forward_euler(integer_order, h, f, x)
+function coefficients = forward_euler(integer_order, xs, x)
   ##
-  ## Use the forward Euler method to compute the derivative of `f` at
-  ## a point `x`.
+  ## Return the coefficients of u(x0), u(x1), ..., u(xn) as a vector.
+  ## Take for example a first order approximation, with,
   ##
-  ## INPUTS:
+  ##   xs = [x0,x1,x2,x3,x4]
   ##
-  ##   * ``integer_order`` - The order of the derivative.
+  ##   f'(x=x1) ~= [f(x2)-f(x1)]/(x2-x1)
   ##
-  ##   * ``h`` - The step size.
+  ## This would return [0, -1/(x2-x1), 2/(x2-x1), 0, 0]. This aids
+  ## the solution of linear systems.
   ##
-  ##   * ``f`` - The function whose derivative we're computing.
-  ##
-  ##   * ``x`` - The point at which to compute the derivative.
-  ##
-
   if (integer_order == 0)
     df = x;
     return;
   end
-    
-  ## We need a few points around `x` to compute the derivative at `x`.
-  ## The number of points depends on the order.
+
   if (even(integer_order))
     offset_b = integer_order / 2;
     offset_f = offset_b;
@@ -31,15 +25,19 @@ function df = forward_euler(integer_order, h, f, x)
     offset_f = offset_b + 1;
   end
 
-  backward_xs = [x-(offset_b*h) : h : x];
-
-  ## We'll always have at least one forward point, so start this vector
-  ## from (x + h) and include `x` itself in the backward points.
-  forward_xs = [x+h : h : x+(offset_f*h)];
+  ## Zero out the coefficients for terms that won't appear. We compute
+  ## where `x` is, and we just computed how far back/forward we need to
+  ## look from `x`, so we just need to make the rest zeros.
+  x_idx = find(xs == x);
+  first_nonzero_idx = x_idx - offset_b;
+  last_nonzero_idx = x_idx + offset_f;
+  leading_zero_count = first_nonzero_idx - 1;
+  leading_zeros = zeros(1, leading_zero_count);
+  trailing_zero_count = length(xs) - last_nonzero_idx;
+  trailing_zeros = zeros(1, trailing_zero_count);
 
-  xs = horzcat(backward_xs, forward_xs);
+  targets = xs(first_nonzero_idx : last_nonzero_idx);
+  cs = divided_difference_coefficients(targets);
 
-  ## Now that we have all of the points that we need in xs, we can use
-  ## the divided_difference function to compute the derivative.
-  df = divided_difference(f, xs);
+  coefficients = horzcat(leading_zeros, cs, trailing_zeros);  
 end
-- 
2.43.2