#!/usr/bin/octave -q

% Jonathan Senning <jonathan.senning@gordon.edu>
% Gordon College
% March 19, 1999
% Revised October 24, 2008
% Revised December 11, 2008 to be compatible with both Octave and Matlab
%
% This program demonstrates the use of Taylor Series to solve the initial
% value problem x'(t) = t*sin(t), x(0) = 1.

% Set desired order of Taylor Series

m = 4;

% Set interval and initial condition

a = 0;
b = 20;
x0 = 1;

% Desired number of subintervals

n = 20;

% Begin the solution process...

h = ( b - a ) / n;
t = a;
x = x0;

% T and X are arrays that will contain each t value and each x value
% computed.  They are not necessary for the final solution, but are used
% to produce a plot of the solution curve.

T = t;
X = x;

% Enter the main loop.  Since the equation we are solving is
% x'(t) = t*sin(t) we will need to evaluate both cos(t) and sin(t)
% repeatedly for each value of t.  We compute them once ahead of time.

for i = 1 : n

  c = cos( t );
  s = sin( t );

  % The following is the pattern for the derivatives of x(t). The
  % subscript in each case is the order of differentiation.  Notice
  % the pattern of negative signs and the uniform increments of the
  % coefficient of the second term.
  %
  %   dx(1) =   t * s;
  %   dx(2) =   t * c + s;
  %   dx(3) = - t * s + 2 * c;
  %   dx(4) = - t * c - 3 * s;
  %   dx(5) =   t * s - 4 * c;
  %   dx(6) =   t * c + 5 * s;
  %   dx(7) = - t * s + 6 * c;
  %
  % In the loop below, s1 and s2 are used to figure the signs of the
  % first and second terms; s0 is used as an intermediate variable.
  % It's kind of a mess, but all this loop does is compute the values
  % of dx in the pattern shown above.

  s0 = 1; s1 = -1; s2 = -1;
  for k = 1 : m
    s0 = -s0;
    s1 =  s0 * s1;
    s2 = -s0 * s2;
    if ( mod(k,2) == 0 )
      dx(k) = s1 * t * c + s2 * (k-1) * s;
    else
      dx(k) = s1 * t * s + s2 * (k-1) * c;
    end
  end

  % Now compute the sum of the first m+1 terms of the Taylor series.
  % Using a nested form of the truncated series saves many operations.

  sum = dx(m);
  for j = m : -1 : 2
    sum = h * sum / j + dx(j-1);
  end
  x = x + h * sum;

  % At this point x contains the value x(t+h) so we need to compute the
  % new t value.  This could be done with either t = t+h or the form below.
  % The from below is preferred to avoid accumulated roundoff error from
  % the repeated addition of h.

  t = a + i * h;

  % Now we can save the current values of t and x(t) for plotting.

  T = [ T, t ];
  X = [ X, x ];
    
end

% Done -- plot to compare with true solution: x(t) = sin(t) - t*cos(t) + 1.

T2 = b * [0:5*n] / (5*n);
Y = sin(T2) .- T2 .* cos(T2) .+ 1;

plot( T, X, 'o', T2, Y, '-' );
legend( 'Taylor Solution', 'Exact Solution', 'location', 'South' );
xlabel( 't' );
ylabel( 'x' );
title( sprintf( ['Taylor Series Solution (Order %d) of x''(t) = t*sin(t), ' ...
                 'x(0) = 1'], m ) );
drawnow();

yesno='';
yesno = input( 'Create PNG file? [y/N] ', 's' );

if ( length( yesno ) > 0 && ( yesno(1) == 'Y' || yesno(1) == 'y' ) )
  pngfile = 'taylor.png';
  print( pngfile, '-dpng' );
  fprintf( 'PNG file %s was created\n', pngfile );
end

% End of File