clear; close all; clc;
%% Part a. Check the function given is a solution
% In order to compare the different numerical methods use to solve
% differential equation, in this first part you should verify the function
% given:
%%
%
% $$Y(t) = 2\,t - \exp(-6t)$$
%
% is a solution of
%%
%
% $$y'(t) = 6\:\exp(6t) * (y(t) - 2t)^2 + 2 ,\quad y(0) = -1$$
%
% Thus, the symbolic approach is used to verify them

syms t y
f1 = 6 * exp(6*t) * (y - 2*t)^2 + 2; % y'(t)
f2 = 2*t - exp(-6*t);             % Y(t)
Yprime = diff(f2);

% Now let check if f1 and f2are related.
% In f1 let us substitute t = 0 and y = -1, while in Yprime let us do it
% with t = 0

val_f1 = vpa(subs(f1,[t, y],[0, -1]));
val_f2 = vpa(subs(Yprime,t,0));

if double(val_f1) == double(val_f1)
fprintf('the function Y(t) is solution of d(y(t)dt)\n')
else
fprintf('the function Y(t) is solution of d(y(t)dt)\n')
end

%% Part b. Applying Euler method
% For three values of h solve the differential equation and compare with
% the exact solution at IVP y(0) = -1.

h = [0.05, 0.1];
t0 = 0;       % Initial time
t1 = 1;       % Final time
y0 = -1;    % Value of y function at t0

solPartB = struct();
for i=1:length(h)
solPartB.(sprintf('H_%d', i)) = euler_h(matlabFunction(f1), ...
   t0, t1, y0, h(i));
end

% This part is just for h = 0.15. The initial condition was set to y(0.01)
% calculated before using h = 0.05

solPartB.('H_3') = [0, -1; euler_h(matlabFunction(f1), ...
   0.1, t1, solPartB.H_1(3,2), 0.15)];

exactSolution = double(vpa(subs(f2,t,1)));

figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');
plot(solPartB.H_1(:,1), solPartB.H_1(:,2),...
'DisplayName','h = 0.05','LineStyle','none',...
'MarkerSize',8,'Marker','v','LineWidth',2);
plot(solPartB.H_2(:,1), solPartB.H_2(:,2),...
'DisplayName','h = 0.1','LineWidth',2,...
'MarkerSize',10,'Marker','o');
plot(solPartB.H_3(:,1), solPartB.H_3(:,2),...
'DisplayName','h = 0.15','LineWidth',2,...
'MarkerSize',8,'Marker','s');
plot(1,exactSolution,'DisplayName','Exact',...
'MarkerFaceColor',[0 0 0],'MarkerSize',15,...
'Marker','pentagram',...
'LineStyle','none',...
'Color',[0 0 0]);
ylabel('y(t)');
xlabel('t');
%ylim(axes1,[-2 10])
box(axes1,'on');
set(axes1,'FontSize',12,'XGrid','on','YGrid','on');
legend1=legend(axes1,'show');
set(legend1,'Location','best')

%% Part c. Applying RK4 method
% For three values of h solve the differential equation and compare with
% the exact solution at IVP y(0) = -1.

solPartC = struct();
for i=1:length(h)

5. A N-by-N diagonal matrix is one of the form d 1 d2 , d N...