#Applying the Runge-Kutta method (This code must be implement on a different file than the main one).

function temp = rk4(func,x,pvi,h)
    K1 = h*func(x,pvi);
    K2 = h*func(x+0.5*h,pvi+0.5*K1);
    K3 = h*func(x+0.5*h,pvi+0.5*K2);
    K4 = h*func(x+h,pvi+K3);
    temp = pvi + (K1 + 2*K2 + 2*K3 + K4)/6;
endfunction

#Main Program.

f  = @(t) (1/16)*((t.^2 + 4).^2);
df = @(t,y) t*sqrt(y);

pvi = 1.0;
h   = 0.1;
Yn  = pvi;

for x = 0:h:10-h
    pvi = rk4(df,x,pvi,h);
    Yn = [Yn pvi];
endfor

fprintf('Time \t Exact Value \t ODE4 Value \t Num. Error\n');

for i=0:10
    fprintf('%d \t %.5f \t %.5f \t %.4g \n',i,f(i),Yn(1+i*10),f(i)-Yn(1+i*10));
endfor