function Tabular = RK2and4Exercise21at6_2(int,x0,h)
% Main Runge-Kutta (RK2 and RK4, compared with Euler and exact)
% Input: Interval [ta,T], 
%       initial condition x0,  and stepsize h
%       number of steps n=T/h
% Output: the approximate solution x, at time T
% Example usage: Tabular = RK2and4Exercise21at6_2([0 0.3],1,0.1)
% the exact solution of z'=t+z, with IC: z(0)=1;
%   is z(t)=-t-1+2exp(t);
% Remark: for the homework assignment, Exercise 23: 
% z'=x-2*z, z(0)=1, the exact solution is z(x)=x/2-1/4+(5/4)*exp(-2*x).
%h = [0.015 0.020 0.025 0.030];
ta = int(1);
finalT = int(2);
n = floor(finalT/h);
    for timestep = h
    [time,u] = euler(int,x0,timestep);    
    [time,y] = RK2(int,x0,timestep);
    [time,x] = RK4(int,x0,timestep);    
%    exact = -time-1+2*exp(time);    %--  Exercise 21
    exact = time/2-1/4+(5/4)*exp(-2*time);  %--  Exercise 23
figure(1)
hold all
plot(time,u,'o','LineWidth',1,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',10);
plot(time,y,'.-.','LineWidth',1,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',10);
plot(time,x,'.','LineWidth',1,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',10);
plot(time,exact ,'-','LineWidth',1,... 
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',10);  

         
            
title('Euler,RK2, RK4 and exact solution, EX 21 at Section 6.2')
legend('forward Euler','RK2 approximation','RK4 approximation','exact','Location','southeastoutside')            
strtimestep = ['step = ',num2str(timestep)];            
text(finalT,x(end),strtimestep,'HorizontalAlignment','left')             
   % pause
    end
Eulersolution = u';    
RK2solution = y';
RK4solution = x';
Tabular = table(time,Eulersolution,RK2solution,RK4solution,exact);
Tabular(1:length(time),:)  ;  
end

function [time,u] = euler(int,x0,h) %%%% --- forward Euler method
ta = int(1);
finalT = int(2);
n = round((finalT - ta)/h);
time=zeros(n,1);
u(1) = x0;
time(1) = ta;
format long 

for i = 1:n
  u(i+1)=u(i)+ h*xdot(time(i),u(i));
  time(i+1) = time(i) + h;
end
end

function [time,y] = RK2(int,x0,h) %%%% --- RK2 method
ta = int(1);
finalT = int(2);
n = round((finalT - ta)/h);
time=zeros(n,1);
y(1) = x0;
time(1) = ta;
format long 

for i = 1:n 
    s1 = h*xdot(time(i),y(i));
    s2 = h*xdot(time(i)+h,y(i)+h*s1);
    y(i+1) = y(i) + (s1 + s2)/2;
    time(i+1) = time(i) + h;
end
end


function [time,x] = RK4(int,x0,h)  %%%% --- RK4 method
ta = int(1);
finalT = int(2);
n = round((finalT - ta)/h);
time=zeros(n,1);
x(1) = x0;
time(1) = ta;
format long 

for i = 1:n
    K1 = h*xdot(time(i),x(i));
    K2 = h*xdot(time(i)+h/2,x(i)+K1/2);
    K3 = h*xdot(time(i)+h/2,x(i)+K2/2);
    K4 = h*xdot(time(i)+h,x(i)+K3);
    x(i+1) = x(i) + (K1 + 2*K2 + 2*K3 + K4)/6;
    time(i+1) = time(i) + h;
end
end

%function z = xdot(t,x)  %the RHS of the ODE in Exercise 21
%z = t+x;
%end
function f = xdot(x,z) % RHS of the ODE in Exercise 23
f = x-2*z;
end
