function hw3
clear all;
close all;
format long g;
options = odeset('stats','on','AbsTol',1e-300,'RelTol',1e-300);

% % Problem1
% [T,Y] = ode45(@prob1,[0 10],0,options);
% plot(T,Y(:,1),'-');
% [p,dummy] = size(Y);
% Y(p)


% % Problem2
% % [T,Y] = ode15s(@prob2,[0 1e10],[1 0 0] ,options);
% % [p,dummy] = size(Y);
% % Y(p,:)
% % plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')


% Problem3
[T,Y] = ode15s(@prob3,[0 200],[1 0 0] ,options);
%plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')
[p,dummy] = size(Y);
Y(p,:);

%
%These are the most accurate answer
%

a = 0.535839131377919;
b = 4.4875872381329e-006;
c = 0.464156381034849;


%
%Below we calculate the relative error at the right end point
%
an = Y(p,1)
bn = Y(p,2)
cn = Y(p,3)

ea = abs(a - an)/abs(a)
eb = abs(b - bn)/abs(b)
ec = abs(c - cn)/abs(c)





function dy = prob1(t,y)
dy = zeros(1,1);% a column vector
dy(1) = 1 * (y(1) - sin(t)) + cos(t);

function dy = prob2(t,y)
dy = zeros(3,1);% a column vector
dy(1) = -1*1e-9*y(1);
dy(2) = 1e-9*y(1) - 1e9*y(2);
dy(3) = 1e9 * y(2);


function dy = prob3(t,y)
dy = zeros(3,1);% a column vector
dy(1) = -0.04*y(1) + 1e4*y(2)*y(3);
dy(2) = 0.04*y(1) - 1e4*y(2)*y(3) - 3*1e7*(y(2))^2;
dy(3) = 3*10^7*(y(2))^2;


 




