% Fourth Order Runge-Kutta method 

clear all;
close all;
clc;

disp('')
disp('Runge-Kutta applied to a scalar initial value problem')
disp('dy/dt(t)=2*y/t + t^2*exp(t), y(1)=0')

kmax = 3;                  % number of experiments
for k=1:kmax 
    Error(k)=0;
    t=1; y=0;            % Initial Conditions
    m=10^k; h(k)=1/m;   %Step Size (h) and number of subintervals (m)
    j=0;
    fprintf('\nApproxs. for grid of %i every %i grid points\n\n ',m,10^(k-1));
    for i=1:m
        k1=h(k)*f(t,y);
        k2=h(k)*f(t+h(k)/2,y+k1/2); 
        k3=h(k)*f(t+h(k)/2,y+k2/2);  
        k4=h(k)*f(t+h(k),y+k3); 
        t= t+h(k); 
        y= y + 1/6*(k1+2*k2+2*k3+k4);
        LErr=abs(y-fexact(t));        %Computing the error of the approx.
        if LErr > Error(k)
            Error(k)=LErr;
        end
        % Producing a shorter vector for plotting
        if mod(i,10^(k-1)) == 0
             j=j+1;
             tv(j)=t;
             yv(j)=y;    
            disp(sprintf('t(%i)= %15.8e y(%i)= %15.8e',...
                 i,t, i,y))        
        end
        
    end
    figure;
    % plot solution
    plot(tv,yv,'o-');
    title(['Approximate solution for: ',num2str(m),' grid points'],'fontsize',13);
end

% print table containing errors for different time steps h:
% table headings:
disp(' ')
disp('       h          Error(Max)');
for k=1:kmax
    disp(sprintf('%13.4e  %13.4e',...
                 h(k), Error(k)));
end

% Following command produce log-log plot of errors and least squares fit
error_loglog(h,Error);