% RK4 for solving IVP
% y'(t) = f(t,y(t)), y(a) = alpha 
% t defined on [a,b]

R =  7.94953123626628;
r = 0.1;
%g = 32.1;
g = 32.1*60*60;
x0=8;
A = @(t,x)(pi*((r + (R-r)*x/8)^2));
f = @(t,x)( (-0.6)*pi*(r^2)*sqrt(2*g)*sqrt(x)/(A(t,x)) );

a = 0;
b = 10;
alpha = 8; % y(a) = alpha

h = 0.1;

N = round((b-a)/h)

t = [a:h:b];
w = zeros(N+1,1);
err =  zeros(N+1,1);
t(1) =a;

w(1) = alpha;

for i=2:N+1 %i=1:N
    k1 = h*f(t(i-1),w(i-1));
    k2 = h*f(t(i-1)+h/2, w(i-1)+ 0.5*k1);
    k3 = h*f(t(i-1)+h/2, w(i-1) + 0.5*k2);
    k4 = h*f(t(i), w(i-1) + k3);
    w(i) = w(i-1) + +(1/6)*(k1 + 2*k2 + 2*k3 + k4);
    t(i) = a+(i-1)*h;
    disp(['it # = ' num2str(i) ' %% ' 'w = ' num2str(w(i)) ] )
end
%format long g
disp('****************************************************')
disp(' t              w')
disp([num2str([t', w, err],8)])
disp('****************************************************')


