% Igor Yanovsky
% Teaching Assistant
% Math 151B
% 04/28/07
%
%#####################################################################
%
% This program solves the system of two ordinary differential equations
% 
%  du1/dt = f1(t,u1,u2),    u1(a) = u10,
%  du2/dt = f2(t,u1,u2),    u2(a) = u20,
% 
% using RUNGE-KUTTA FOURTH-ORDER method.
% Functions defining the system of ODEs, as well as their exact solution, 
% are defined as separate Matlab functions.
%
%#####################################################################

clear all

format long

a =  0.0;
b =  1.0;
u10 = -1.0;
u20 =  0.0;

h = 0.1;
N = (b - a)/h;

t = a;
w1(1) = u10;
w2(1) = u20;
u1e(1) = u10;   % exact solution
u2e(1) = u20;
err1(1) = 0;
err2(1) = 0;


for(i=1:N)
    
    % advances the system of two ODEs using RK-4:
    k11 = h*f1( t        , w1(i)          , w2(i)           );
    k12 = h*f2( t        , w1(i)          , w2(i)           );
    
    k21 = h*f1( t+(h/2.0), w1(i)+(k11/2.0), w2(i)+(k12/2.0) );
    k22 = h*f2( t+(h/2.0), w1(i)+(k11/2.0), w2(i)+(k12/2.0) );
    
    k31 = h*f1( t+(h/2.0), w1(i)+(k21/2.0), w2(i)+(k22/2.0) );
    k32 = h*f2( t+(h/2.0), w1(i)+(k21/2.0), w2(i)+(k22/2.0) );
    
    k41 = h*f1( t+h      , w1(i)+k31      , w2(i)+k32       );
    k42 = h*f2( t+h      , w1(i)+k31      , w2(i)+k32       );
    
    w1(i+1) = w1(i) + (1.0/6.0)*( k11 + 2.0*k21 + 2.0*k31 + k41 );
    w2(i+1) = w2(i) + (1.0/6.0)*( k12 + 2.0*k22 + 2.0*k32 + k42 );
    
    t = t + h;
 
    u1e(i+1) = u1Exact( t );
    u2e(i+1) = u2Exact( t );
    
    err1(i+1) = abs( w1(i+1)-u1e(i+1) );
    err2(i+1) = abs( w2(i+1)-u2e(i+1) );
end


fprintf(1,'N = %3d',N);  fprintf(1,', h = %4.7e',h);  fprintf(1,', t = %4.2f',t);  fprintf(1,',\n');



fprintf(1,'w1 = %5.10e', w1(N+1)    );  fprintf(1,',\n');
fprintf(1,'w2 = %5.10e', w2(N+1)    );  fprintf(1,',\n');
fprintf(1,'u1 = %5.10e', u1e(N+1));  fprintf(1,',\n');
fprintf(1,'u2 = %5.10e', u2e(N+1));  fprintf(1,',\n');
fprintf(1,'err1 = %5.10e', err1(N+1) );  fprintf(1, '\n');
fprintf(1,'err2 = %5.10e', err2(N+1) );  fprintf(1, '\n');


clear all