#!/usr/bin/octave -q

% Jonathan Senning <jonathan.senning@gordon.edu>
% Gordon College
% March 19, 1999
% Revised April 27, 2006 to update graphing commands
%
% This program demonstrates the use of Euler's Method to solve the initial
% value problem x'(t) = t*sin(t), x(0) = 1.

% 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.

for i = 1 : n

  x = x + h * t * sin(t);

  % 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( 'Euler Solution', 'Exact Solution', 'location', 'South' );
xlabel( 't' );
ylabel( 'x' );
title( sprintf( "Euler Solution (n = %d) of x'(t) = t*sin(t), x(0) = 1", n ) );
drawnow();

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

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

% End of File