%
% Written by:
% -- 
% John L. Weatherwax                2005-08-14
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

clear functions; 
%clear; 
close all; 
clc; 

randn('seed',0); rand('seed',0);

t_final = 1.0;


% Part (a) using the Chandrasekhar type algorithm: 
%
[t_a3,y_a3] = ode45( @sect_5_3_ode_fn_part_a, [0, t_final], [0;sqrt(10);0;0] );

% Part (a) using the method from Problem 5.1:
%
[t_a1,y_a1] = ode45( @sect_5_1_ode_fn_part_a, [0, t_final], [0;0;0] );
K_C_a = [ y_a1(:,2) y_a1(:,3) ];

% Compare the two components of K_C computed two different ways:
figure;
plot( t_a3, y_a3(:,3), '-k' ); hold on;
plot( t_a1, K_C_a(:,1), '-r' ); title('K_C_a(1)');

figure;
plot( t_a3, y_a3(:,4), '-k' ); hold on;
plot( t_a1, K_C_a(:,2), '-r' ); title('K_C_a(2)');

pause;
close all; 


% Part (b) using the Chandrasekhar type algorithm: 
%
[t_b3,y_b3] = ode45( @sect_5_3_ode_fn_part_b, [0, t_final], [0;sqrt(10);0;0] );

% Part (b) using the method from Problem 5.1:
%
[t_b1,y_b1] = ode45( @sect_5_1_ode_fn_part_b, [0, t_final], [0;0;0] );
K_C_b = [ y_a1(:,1) y_a1(:,2) ];

% Compare the two components of K_C computed two different ways:
figure;
plot( t_b3, y_b3(:,3), '-k' ); hold on;
plot( t_b1, K_C_b(:,1), '-r' ); title('K_C_b(1)');

figure;
plot( t_b3, y_b3(:,4), '-k' ); hold on;
plot( t_b1, K_C_b(:,2), '-r' ); title('K_C_b(2)');

pause; 
close all; 


% Part (c) using the Chandrasekhar type algorithm: 
%
[t_c3,y_c3] = ode45( @sect_5_3_ode_fn_part_c, [0, t_final], [0;sqrt(10);0;0;0] );

% Part (c) using the method from Problem 5.1:
%
[t_c1,y_c1] = ode45( @sect_5_1_ode_fn_part_c, [0, t_final], [0;0;0] );
K_C_c = [ y_c1(:,1) y_c1(:,2) y_c1(:,3) ];

% Compare the two components of K_C computed two different ways:
figure;
plot( t_c3, y_c3(:,3), '-k' ); hold on;
plot( t_c1, K_C_c(:,1), '-r' ); title('K_C_c(1)');

figure;
plot( t_c3, y_c3(:,4), '-k' ); hold on;
plot( t_c1, K_C_c(:,2), '-r' ); title('K_C_c(2)');

figure;
plot( t_c3, y_c3(:,5), '-k' ); hold on;
plot( t_c1, K_C_c(:,3), '-r' ); title('K_C_c(3)');